Flow induced by a line sink near a vertical wall in a fluid with a free surface, Part II: finite depth

Flow caused by a line sink near a vertical wall in an otherwise stagnant fluid with a free surface is studied. A linear solution for small flow rates is obtained and a numerical method based on fundamental singularities techniques is applied to the full nonlinear problem. The sink is located at an arbitrary location away from all boundaries and the fluid is of finite depth. Steady solutions are presented for various flow rates and sink location. It is shown that the numerical results and linear solutions are in good agreement for small flow rates. The results suggest that steady nonlinear solutions are limited to flow rates below some critical value. Some interesting surface shapes are obtained depending on the location of the sink.

of waves upstream of the sink, since the flow slows to zero in the far field [6]. In the case of finite depth considered herein, waves are indeed possible. Important work using exponential asymptotics, such as that of Chapman and Vanden-Broeck [7] and more recently Lustri et al. [8], considered the existence of steady waves due to a line source in a fluid of finite depth, a problem which is mathematically identical to the flow into a line sink. The latter showed that waves of exponentially small order exist for the problem of steady flow from a line source as the flow rate approaches zero. These waves are exponentially small, and so will not appear in any regular algebraic power series and are too small to appear in the numerical solutions provided herein.
The behaviour of such flows is significant in managing water quality in different water storages and cooling ponds [9] and led to the publication of many papers discussing withdrawal through a line sink beneath a free surface and related problems. The model used here is also useful as an idealisation of withdrawal from a two-layer, density stratified fluid if the interface between the layers is thin and one of the layers is almost stagnant, such as in lakes or reservoirs.
Two types of steady solution exist for a single-layer flow in the absence of surface tension. The first type (see for example [6,[10][11][12][13]), are those for which the free surface pulls down into a cusp directly above the sink. Hocking [14,15] showed that these solutions are the point of transition to two-layer flows such as those described in [16], in the case of a fluid with a two-layer stratification. If the lower layer is unbounded these solutions occur at a unique flow rate for a given geometry, while if the lower layer is bounded below, Vanden-Broeck and Keller [17] showed that cusp solutions exist for a range of different flow rates with F > 1 and that there is also a set of unique solutions for F < 1.
The other type of single-layer flow is one in which a stagnation point forms on the free surface above the line sink [1,2,11], and similar solutions are described for a point sink in Hocking et al. [18]. In both of these cases, steady solutions exist for all values of flow rate beneath some critical value. Forbes and Hocking [1] found the interesting result that when surface tension is included, there is a bifurcation in the parameter space leading to multiple different solutions for the same geometric configuration.
Experimental work [19][20][21][22] showed considerable variation in the values of the critical flow rate at which the transition from single-layer to two-layer flow occurs over a range of different geometries. A theoretical investigation of the flow is important in interpreting results obtained in experiments and in the field.
Here, we consider a single-layer flow in which a line sink is situated in an arbitrary location near a vertical wall situated at x = 0. The problem is characterised by several nondimensional quantities, the first of which is the Froude number, defined to be where g is gravity, U is the upstream velocity, m is the flux from the region into the sink, and H B is the depth of the layer beneath the upstream level of the interface. Note that the velocity U = m/H B , where m is the flux into the line sink. Other important quantities are the sink location, i.e. (x s , y s ), and surface tension which is defined as where T is the surface tension parameter. Almost all previous work has focussed on a sink located on the wall, and so this work considers the impact of a nearby wall on the critical parameters of the flow.

The governing equations
Consider the steady two-dimensional, irrotational flow of an inviscid, incompressible fluid toward a line sink situated at an arbitrary location beneath a free surface y = η(x) as shown in Fig. 1.
The fluid is assumed to have a free surface and the shape of it will be examined both with and without surface tension, T . The problem can be formulated in terms of a velocity potential φ(x, y), where the horizontal velocity component is u = φ x and the vertical velocity component is v = φ y , where x is the horizontal coordinate and y is the vertical coordinate. Nondimensionalising the velocity and the length with respect to U = m/H B and H B , respectively, the problem is to solve subject to the dynamic condition obtained from setting pressure to the atmospheric value on the free surface via the Bernoulli equation, i.e.
where u = (u, v) is the velocity vector. Note that in this formulation the surface η(x) → 1 as x → ∞. There are also kinematic conditions that ensure there is no flow through the free surface or the wall in steady flow, where the subscript denotes partial differentiation, y = b(x) is the equation of the bottom surface and y = η(x) is the equation for the elevation of the free surface. In this work we let b(x) = 0 for all x. In this formulation, information is lost as the sink height, y s approaches 1 and so it is convenient to define a modified Froude number based on the depth of the sink rather than the depth of the fluid, taking the flux-defined form of F from (7). In the limit as we approach the line sink at (x, y) = (x s , y s ) the velocity potential should take the form which corresponds to a sink of unit strength, i.e. there is one unit of fluid entering per unit time.

The linearised problem
It is possible to find an approximate solution to this problem if the flow rate is assumed to be quite small (F << 1) and hence the disturbance of the surface is small. Much can be learned from such solutions and they also provide a strong verification of the numerical solutions computed in the next section. An expansion of the dynamic condition about the flow along the top of the duct, if the top of the duct is assumed to be flat (the rigid-lid assumption) is used to approximate the free surface.
To compute such a solution, we let and assume the free surface disturbance is small so that on the surface all of the terms are small. All quantities are evaluated on y = 1. By simplifying and rearranging the dynamic and kinematic condition, we obtain the equations for the leading order problem as and

The flow in a duct
The leading order problem (9)(10)(11) represents the flow in a duct as shown in Fig. 2 with the surface at AI being horizontal at y = 1, i.e. Y 0 = 1. Conformal mapping techniques can be used to map the flow in the upper half-plane ζ = ξ + iχ to the real duct in the physical z = x + iy plane in order to satisfy Eqs. (9)(10)(11). The appropriate mapping from the z-plane to the ζ -plane is given by We choose to let ζ s = (ξ s , χ s ) correspond to the line sink in the ζ -plane, so that the boundary AI maps to the line z = x + i, 0 ≤ x < ∞, since the free surface lies on −∞ < χ < −1. The complex potential, f 0 (ζ ) for the flow due to a line sink at ζ = ζ s , with no flow through the real line in the ζ -plane, where the location of the sink in the z-plane is such that ζ s = cosh π z s , is where ζ denotes the complex conjugate, so that the complex potential in the z-plane is log(cosh π z − cosh π z s ) + log(cosh π z − cosh π z s ) .
representing the velocity potential and the streamfunction, respectively, are and where A = cosh 2 π x cos 2 π y − sinh 2 π x sin 2 π y + 2 cosh π x cosh π x s cos π y s cos π y + cosh 2 π x s cos 2 π y s + sinh 2 π x s sin 2 π y s , and B = 2 sinh π x sin π y cosh π x s cos π y s + 2 cosh π x cos π y sinh π x sin π y.
Streamlines, ψ 0 (x, y), from (19) and contours of the velocity potential, φ 0 (x, y), from (18) for several different sink locations are shown in Figs. 3 and 4. As the sink moves away from x = 0 a region of very slow flow occurs close to x = 0, which can be seen as the large separation of streamlines in Fig. 3 or a long distance between equipotentials in Fig. 4. This region grows in size as x s increases. Computing all quantities along y = 1 in (13), the approximation at O(F 2 ) yields a second-order ordinary differential equation in Y 1 (x), where Y 1 (x) is the first-order perturbation to y = η(x), and φ 0,x (x, 1) can be determined from (19). The variation of parameters method can be used, as in [1] for flow due to a line sink in water of finite depth, and the approximate solution for η(x) to order O(F 2 ) can be computed as where G(x) = φ 2 0x (x, 1). These integrals are not easily solved exactly, so we have evaluated them numerically using octave routine quad. Let us first consider the case with no surface tension.
Typical free surfaces are shown in Fig. 5, for several different sink heights, y s = 0.2, 0.6, 0.8, with x s = 2 and F = 0.1. There appear to be 3 different free surface shape types depending on the location of the sink. When y s is located far away from the surface, y s = 0.2, the surface is flat downstream and upstream and decreases monotonically as x increases. As y s moves up from the bottom, y s = 0.6, a trough is formed on the upstream side of the sink location while on the wall side it remains flat. As y s moves up further and gets near the surface, y s = 0.8, a peak forms on the free surface as the surface rises back up to y = 1 + 1 2 F 2 , the stagnation level from Eq. (3), almost directly above the sink, and a second trough forms on the wall side along with the one on the upstream side. Thus, Fig. 5c shows that a second stagnation point has formed on the free surface, and not in the other two cases. The shape of the free surface depends on the location of this stagnation point in the flow domain. Note that in the absence of surface tension there is always a stagnation point at x = 0, y = 1 + 1 2 F 2 , but if there is a second stagnation point on the upper surface, then a trough and a peak will form as described above. The location and circumstances of this second surface stagnation point can be determined from the mapping. Note that there is always a stagnation point located at the base of the wall at (0, 0). These shapes can be determined from the properties of the mapping and so we have the following important result. The approximate solution predicts a third stagnation point in the domain as follows; ⇒ stagnation point on the bottom.
⇒ stagnation point on the side wall. Under the conformal transformation (15), as shown in Fig. 2, if ξ lies on −∞ < ξ < −1 then it maps to the line AI , the free surface, if ξ lies on −1 < ξ < 0 it maps to the line I C, the vertical wall, and if ξ > 0 it maps to the line C J , the channel base. Therefore, if Re{cosh π z s } < −1, then z * = z s lies along z = x + i and hence will be on the free surface. If −1 < Re{cosh π z s } < 0 then z * lies along z = iy and hence will be on the vertical side wall at x = 0. If Re{cosh π z s } > 1 then z * lies along z = x + 0i and hence will be on the bottom.
These regions are indicated in Fig. 6 in the z-plane, with sink locations (x s , y s ) above the upper curve being all of the possible locations for which the extra free-surface stagnation point forms. The lower curve in Fig. 6 shows the upper bound on the region in which the sink produces a stagnation point on the bottom of the channel. If the secondary stagnation point occurs on the free surface, it indicates a flow convergence at this point, with fluid coming from both upstream and downstream of the sink location before being drawn down and out.

Solution using the fundamental singularities method
To consider the full nonlinear steady flow problem we need to implement a numerical scheme. We have chosen a fundamental singularities method such as that first used by Scullen and Tuck [23] to compute surface waves, and subsequently by Stokes et al. [4] and Mansoor et al. [5] for a line sink, a method also implemented by Stokes et al. for a point sink [24]. The main idea is that we assume a set of source or sink points in some positions (X j , Y j ) just above the free surface, and compute the velocity potential function φ(x, y), for the existing surface shape. The points are assumed to be symmetrically distributed about the x and y-axes, and all sources and sinks are likewise assumed to be symmetrically distributed in terms of strength.
We take a discrete form (x k , η k ), k = 1, 2, . . . , N of the free surface. A function φ that satisfies Laplace's equation (2) and the condition of no flow through the boundaries (4) and (6) is where φ 0 is the solution for flow in a duct given in (18), γ k are strengths of the sinks situated just above the free surface, and are the distances between the kth source or sink point and the point (x, y), which lies in the fluid. The sink strengths γ k and surface heights η k are unknown and will be determined once the value of φ is known at the N distinct surface points, (x j , η j ), j = 1, 2, . . . , N .
The locations (X j , Y j ) need to be quite close to the surface (but outside of the flow domain) so that a wellconditioned system of equations is derived. It was found that choosing the locations at a distance of one-third the local grid separation, and perpendicular to the free surface as in Scullen and Tuck [23] or in Stokes et al. [4], provided accurate solutions. Trial and error indicated that the solutions were not particularly sensitive to this choice.
The computational domain was truncated at x = L. The choice (25) satisfies conditions (2), (4) and (6) using the method of images which is why there are 4 log terms present. If the number of discrete points used to represent the surface is N , the η k values and the coefficients γ k of the sink term can be used as unknowns at the fixed values of x k , k = 1, 2, . . . , N . The equations to be satisfied are the steady components of Bernoulli's equation (3), and the kinematic condition (5), both to be satisfied at each of the N points. This gives 2N equations in 2N unknowns, and an iterative Newton's method was implemented to solve for the unknowns. On most occasions, starting guesses for the sink strengths were set to zero and the surface height was set at y = 1. This method was programmed in Octave and also Fortran to enable calculation with very large values of N . Convergence was rapid (usually within 4 or 5 iterations) except as limiting values were approached. Solutions converged to graphical accuracy in most cases with N = 300, and most calculations were done with N = 500.
As the value of F increased, the numerical method began to exhibit grid-scale oscillations and failed to converge in the case with zero surface tension β = 0. This not surprising as similar behaviour has been seen in most previous, related works [1,2,11]. This was found to occur for all sink locations and was the maximum value of Froude number, F max , at which we could compute solutions. As the maximum F solutions were approached, higher values of N = 1000 were used to verify the limiting values, but the behaviour was the same. In cases with nonzero surface tension, these oscillations were not observed, but eventually the numerical scheme failed to converge in these cases also.

Results
After computing a large number of solutions, it became clear that all of the results obtained for this problem depend on three factors, the Froude number F, the sink location (x s , y s ) and the surface tension β. We begin with the results of what was observed in the zero surface tension case β = 0. Another comparison of the linear and nonlinear solutions is given in Fig. 8, but with the sink located at (x s , y s ) = (2, 0.8). It can be seen that the two results are again in good agreement, although nonlinear effects cause a slight deepening of the surface in the bigger trough. An investigation showed that the predictions of the shape of the surface relative to sink location from the linear solution, are closely followed by the nonlinear numerical solutions. The second stagnation point on the surface for the case (x s , y s ) = (2, 0.8) can be seen as the surface rises back to the stagnation level y = 1 + 1 2 F 2 . No other free-surface stagnation points exist. The very slow flow region between the sink and the wall can be seen as the very flat level of the surface.
As the Froude number increases, the depth of the main trough increases as expected and we also observe a relatively short narrow peak on the wall side of the trough. The free surface is flat in the far field with no waves observed. A nearly stagnant zone is again evident near the wall, but a second stagnation point is now visible on the surface between the sink and the wall, where the surface again rises to the stagnation height y = 1 + 1 2 F 2 when β = 0.
The effect on the surface shape of a small amount of surface tension β = 0 can be seen in Fig. 9, which shows a comparison between the numerical and linear solutions with the sink at (x s , y s ) = (2, 0.8) and β=0.001. The added surface tension term has made very little difference to the shape of the surface, but even this small amount means that the nonlinear solutions can be calculated to a significantly larger value of Froude number, F ≈ 0.2 (c.f. F ≈ 0.13 with β = 0). Even at this larger value of F there is only a slight difference between the two solutions, except very close to the limiting case of the nonlinear solution, when the dip at x ≈ 2.3 narrows and deepens dramatically. The nonlinear solution has a narrowing and great deepening in the trough above the sink, as the maximum F max ≈ 0.2 solution is approached. This rapid deepening is only observed near the maximum flow rate when surface tension is included.  These results demonstrate clearly that the linear solution represents the shape of the free surface very well even for quite large values of flow rate. At the same time it verifies that the numerical scheme is working correctly.
In all cases the linear solution provides a very good depiction of the free surface shape but it can not be used to determine the maximum Froude number for each case. The maximum solutions including surface tension were easier to obtain as the grid-scale oscillations did not occur due to the regularising effect of surface tension. As surface tension β increased, the maximum computable solution occurred at larger values of F. Figures 11 and 12 show the maximum Froude number computed for different y s and β and x s = 1. When β = 0 the highest value F max is approximately F max ≈ 0.246 when the sink is close to the bottom and it remains at around this value until y s ≈ 0.5 (see Fig. 11). However, F Smax , the maximum Froude number dependent on the depth of the sink instead of the fluid depth, is approximately F Smax ≈ 0.69 if the sink is close to the bottom and steadily increases to a maximum at y s ≈ 0.9, where it is F Smax ≈ 1.391 before decreasing again (see Fig. 12). This parameter F S is more representative of the effect of the sink on the surface, and its values are consistent with the Froude number defined in the first of these two articles [5].
At any fixed sink location the maximum values of F and F S increase as surface tension β increases, as shown in Figs. 11 and 12. This rise is very rapid, even for small β. For the larger value of surface tension β = 0.1, the maximum is getting larger up to F max ≈ 0.9 when the sink is located very close to the bottom and F max remains approximately constant until y s = 0.7 then drops to zero as y s → 1 (as the sink approaches zero depth). However, F S rises from F Smax ≈ 0.9 when the sink is close to the bottom to a maximum around F Smax ≈ 7.6 when y s ≈ 0.8 before dropping slightly as the sink approaches the surface. If the value of the surface tension β = 0.2 the maximum is much higher again. Figures 11 and 12 show that the maximum is F Smax ≈ 1 where the sink is close to the bottom and F Smax ≈ 8 if the sink is close to the surface. At higher values of surface tension the maximum value of F S remains quite similar until y s ≈ 0.75 before dipping downward. In Fig. 12, this is reflected by the peak in F Smax occurring at y s ≈ 0.8. This is to be expected since the surface tension stabilises the surface, enabling steady flows to occur at higher flow rates. Forbes and Hocking [1] showed that this maximum corresponded to a bifurcation in solution space and nonunique solutions.
In Figs. 13 and 14 the surface tension value is increased to β = 0.01 and β = 0.1. It is very clear that the effects are mostly near the peak, which is getting lower, while the trough is getting shallower. Note that the secondary surface stagnation point has disappeared.
A comparison of solutions with different surface tension values is shown in Fig. 15. The three solutions were computed using the same Froude number F = 0.14 and the sink is close to the surface (x s , y s ) = (1, 0.8), with different surface tension values β = 0.001, β = 0.01 and β = 0.1, respectively. As observed in Fig. 15, the three solutions are different. If β = 0.001 there is a slightly shorter peak on the wall side and a deeper trough. However, with an increase in the surface tension value, the peak disappears and the trough gets shallower. In the case with the sink at y s = 0.8 and x s = 1 and surface tension as small as β = 0.001 the limiting fundamental singularities method solution occurs when F max ≈ 0.68. The horizontal location of the sink has a big influence on the maximal flow rate at which steady solutions can be found until y s gets close to the surface. As x s increases it would be expected that the shape should approach that of an isolated sink (with no wall present). However, even for quite large values of x s the surface shape near the sink is quite asymmetric. The values of F Smax approach the values of the equivalent Froude number in Part I of these two articles [5] as y s → 1.
These computations were made for a number of different configurations and the following trends were observed. Raising the sink closer to the surface decreased the maximum flow rate at which a steady solution could be obtained (see Figs. 16 and 17). These figures also show the dramatic rise in F max when surface tension is increased. This is most dramatic once the sink is very close to the free surface such as 0.7 < y s < 1. Decreasing the value of y s leads to an increasing maximum flow rate. That this general trend should be apparent is not surprising since the sink needs to "pump" harder to get the same distortion to the surface.
The effects of moving the sink away from the wall, shown in Figs. 16 and 17, turn out to be significant with regard to the maximum rate. At constant β and y s , the maximum Froude number increases slightly as x s increases. When close to the wall, the maximum F is lower because the flow is constrained by the wall and as a consequence the sink appears stronger. In fact, as the sink approaches the wall it approaches double the strength (by the method of images). As x s increases the value of F max levels off toward values obtained in the case of an isolated sink in a fluid of finite depth [16].

Concluding remarks
The flow generated by a line sink at an arbitrary location in a channel of finite depth has been considered, and linear and numerical solutions have been computed for the shape of the free surface. Our main focus was the effect of the wall when compared to flows where the sink is actually in the end wall or completely isolated. As the sink moved away from the wall it was found that the solutions obtained moved between these two limiting cases. However, at moderate values of sink separation, x s , it was found that a large, almost stagnant region formed between the outlet and the wall (see Figs. 3, 4), a finding that has important implications for water management in reservoirs, for example. Most of the behaviour is demonstrated in the approximate solution of Sect. 3.2, but a full nonlinear, numerical solution was required to obtain the maximal steady solutions.
This work also showed that there are three different types of surface shape; one with dips and a single stagnation point on the wall, one with a single dip on the wall side of the sink location and hence two stagnation points, and one with dips on both sides of the sink, all dependent on the physical location of the sink.
Lustri et al. [8] computed waves of exponential order for the related problem of a source in a channel of finite depth. Neither the linear nor the nonlinear solutions would pick up these very small waves, but it is possible that their existence may be one cause of the failure of the numerical scheme.
As in earlier work the addition of surface tension allowed for a more thorough investigation into the characteristics of the solution. Even a small amount of surface tension resulted in a solution for larger Froude number, and it also had a big effect on the shape of the free surface.
For nonzero surface tension β = 0, we find that the maximum flow rate, F max , increases to quite a large value, around three times higher with β = 0.05 and 5 times higher with β = 0.2, and as x s increases the maximum Froude number increases. The Froude number based on sink depth F S increases as the sink is situated higher in the fluid, but does seem to reach a local maximum when the sink elevation is 80-85 % of the fluid height.