Shallow free-surface Stokes flow around a corner

The steady lateral spreading of a free-surface viscous flow down an inclined plane around a vertex from which the channel width increases linearly with downstream distance is investigated analytically, numerically and experimentally. From the vertex the channel wall opens by an angle α to the downslope direction and the viscous fluid spreads laterally along it before detaching. The motion is modelled using lubrication theory and the distance at which the flow detaches is computed as a function of α using analytical and numerical methods. Far downslope after detachment, it is shown that the motion is accurately modelled in terms of a similarity solution. Moreover, the detachment point is well approximated by a simple expression for a broad range of opening angles. The results are corroborated through a series of laboratory experiments and the implication for the design of barriers to divert lava flows are discussed. This article is part of the theme issue ‘Stokes at 200 (Part 1)’.


Introduction
Viscous gravity currents are abundant in nature and industry and their modelling began with Sir George Stokes' equations for creeping flow. His ideas in this area have been applied to gravitationally driven viscous flows in diverse fields including glass manufacture, lava flows, coating and printing processes, food manufacture and many biological settings [1]. These flows often encounter fixed obstructions, for example in thin film flow around nuts within aeroengine chambers or when a  Figure 1. Schematic for gravitationally driven viscous flow down an inclined plane at an angle β to the horizontal. The channel expands with opening angle α to the downslope direction, which we take to be the X axis, while the Y axis is in the cross-slope direction. (Online version in colour.) lava flow interacts with a barrier [2,3]. Detailed knowledge of the deformation of the free surface is essential for understanding and informing engineering decisions and design. The interaction between gravity-driven flows and topographic mounds has been well studied and controls on when the mound is surmounted by the flow have been determined [4][5][6][7][8][9]. There have also been theoretical, numerical and experimental analyses of the flow behaviour upstream of surface piercing cylinders of various cross-sections [2,[10][11][12]. These works have shown that there may be a dry zone (in which there is no fluid) downstream of the obstruction and there have been investigations into the location and extent of the dry zone in the regime that capillary forces play a dominant role [11]. However, at larger length scales, for which the effects of surface tension are negligible, the dependence of the dry zone shape on the obstruction cross-section is not yet well understood.
The present work addresses aspects of this problem and our results have important implications for the optimal design of lava flow barriers. While the previous upstream analyses provide valuable estimates of the barrier height that is required to prevent overtopping by the oncoming flow, it is also important to quantify the controls on the downstream dry zone. The downstream behaviour is challenging to model because of the complex interaction between the spreading fluid and the back of the obstruction, which may lead to the fluid detaching from the barrier to leave a 'dry' (fluid-free) zone. To make progress in studying the downstream interaction, we consider the more straightforward geometry of a semi-infinite channel whose wall opens at a vertex downstream (figure 1). The flow detaches from the angled wall at some distance downstream and a major aim of this work is to quantify how the detachment distance along the wall, d, depends on the wall opening angle, α and the depth of the oncoming flow. It is noteworthy that the viscously dominated dynamics analysed in this study are very different from situations in which the flow has significant inertia and detaches at the vertex.
Since a constant flux of fluid is supplied from upstream, the only length scale in the problem is the oncoming flow depth, H ∞ . This enables the derivation of a parameter-free governing equation in §2, using a lubrication approximation. We introduce the numerical method in §3 and present contour plots of the steady flow depth. The dimensionless wall detachment location depends only on the opening angle, α. We show that the contact line where the flow depth vanishes is relatively insensitive to the opening angle and is well described by a simple similarity solution, with a single virtual origin for any opening angle. The results demonstrate that away from the vertex, the influence of the wall and the downslope component of gravitational slumping are negligible and the flow is governed by a balance between the component of gravity acting down the slope and the cross-slope slumping ( § §4 and 6). In addition, we construct an exact solution for the special case of α = π , which corresponds to the abrupt end of the channel wall at x = y = 0 ( §5). This solution is used to validate the numerical results and to demonstrate the utility of the similarity solution. New experimental results are presented in §7 and these show good agreement with our theoretical and numerical predictions for the detachment distance. We conclude with a brief discussion of the implications of our results in the context of lava flow barrier design and infer that the design of the downstream side of barriers has little influence on the extent of the fluid-free zone that is safely protected.

Formulation
We consider flow down a channel inclined at an angle β to the horizontal (figure 1) and focus on the steady state that develops at long times after transient effects associated with the passage of the front have decayed. The channel opens at the origin in (X, Y) coordinates so that the shape of its wall is given by Y = 0 for X < 0 and Y = X tan α for X > 0 (figure 1). Fluid is supplied with constant flux, Q per unit width, from a line source far upstream, which is semi-infinite in Y (−∞ < Y < 0). The flow thickness far upstream is a constant that may be obtained by balancing the downslope component of gravity with the viscous stresses [13,14] where μ, ρ and g are the dynamic viscosity, the density of the fluid and gravitational acceleration, respectively. We use the lubrication approximation to determine the hydrostatic pressure and velocities in the fluid [15]. Local mass conservation is applied to the velocities to obtain [14] ∂H In this expression we have assumed that the effects of surface tension are negligible; for the experiments described below the capillary length is approximately 2 mm [12], which is much smaller than the streamwise length scale. On adopting dimensionless variables h = H/H ∞ and (x, y) = [X/(H ∞ cot β), Y/(H ∞ cot β)], we find that the steady governing dimensionless equation is while the dimensionless flux in the x and y directions is We next describe the boundary conditions for the steady flow. The line source is far upstream of the wall vertex and the line source is semi-infinite. Thus h → 1 as x → −∞ and h → 1 as y → −∞. There is no-flux, q · n = 0, into the angled wall, where the normal direction to the angled wall is n = (sin α, − cos α). Thus the boundary condition on y = x tan α, x > 0, may be written as Downstream of the point where the flow detaches from the wall, there is a contact line, y = y c (x), where the current depth first vanishes: h(x, y c (x)) = 0. We denote the dimensionless distance along the wall at which the flow detaches by d. This nondimensionalized distance depends only on the opening angle α and we write d = f (α). One of the main aims of this work is to determine this dependence. We note that f is a decreasing function of α and f (α) → ∞ as α → 0 because the flow does not detach in this limit.
Our model neglects the no-slip boundary condition at the channel wall in accordance with the leading order lubrication model. We have also neglected surface tension, which may become very significant near the contact line (see, for example, [16]). A detailed local treatment of these effects would be needed to complete the solution but they do not affect bulk behaviour and hence can be neglected for our purposes. For a detailed study of the flow near the wall associated with no-slip and matching with the bulk flow, see [17].

Numerical method
We use a finite-element method to integrate the steady governing equation (2.3) numerically. The routine is part of Matlab's Partial Differential Equation Toolbox TM , which uses adaptive mesh generation. The steady state is found iteratively; we take an initial guess to be h = 1 everywhere and iterate until a converged solution is found. A similar approach has previously been used for flow over a mound and flow past a cylinder [9,12]. We solve the system with the following boundary conditions. There is no normal flux (q · n = 0) into the wall along y = 0 for x < 0 and y = x tan α for x > 0. There is also no flux across the domain boundary at y = −c. Constant flux (q · n = 1) is supplied at the upstream boundary, x = −a, −c < y < 0. Finally, we apply a free-flux condition (∂h/∂x = 0) at the downstream boundary at x = b, −c < y < b tan α. The domain size is determined by the parameters a, b and c. The system is solved with an initial choice of a, b and c and then they are increased and the system solved again. This process is repeated until increasing the domain size further has a negligible influence on the solution. In the case α = π/4, for example, the domain size is given by a = 5, b = 20 and c = 15. Note that the contour plots presented in this paper include only the region of interest rather than the entire domain of the numerical method. Typically, the mesh for the finite-element method contains approximately 100 000 triangles, which provides sufficient resolution. Contour plots of the steady flow thickness are shown in figure 2 for two values of the angle, α.
The nonlinear diffusive terms on the right-hand side of (2.3) become negligible when h is very small. The original numerical scheme was not effective in this regime. It was necessary to add a small flux out of the open wall to coat the region in which h = 0 with a thin film. The modified boundary condition along y = x tan α is q · n = where the parameter was selected to be as small as possible while enabling the code to run efficiently and typical values were 10 −7 . The influence of is predominantly seen at the steep edges of the current, which are slightly smoothed (for example, see figure 3a). Beyond the steep edge of the current, where we anticipate that the flow depth vanishes, the calculated flow depth is at most h = 10 −2 . We use this value (rather than h = 0) to define the edge of the current and obtain the detachment point and the location of the contact line [y = y c (x)].

Far-field similarity solution (x 1)
Sufficiently far downslope from the wall detachment point, (x 1), the governing equation (2.3) may be approximated by [9,14,18] ∂h 3 ∂x = ∂ ∂y h 3 ∂h ∂y (4.1) and the downstream flow is unlikely to be effected by the details near the vertex. The upstream boundary condition at x = 0 is given by h = 1 in y < 0 and h = 0 in y > 0. Also h → 1 as y → −∞ because the upstream line source is semi-infinite. The system is self-similar with scaling y ∼ x 1/2 and the omitted streamwise diffusive term, ∂(h 3 ∂h/∂x)/∂x, relative to the advective term ∂h 3 /∂x is therefore of order 1/x and negligible in the far-field (x 1). We transform equation (4.1) using the similarity variables, η = y/x 1/2 and h(x, y) = χ (η), to with boundary conditions χ (η 0 ) = 0 and χ → 1 as η → −∞. To solve this system numerically, we shoot from χ = 0 and iterate to determine η 0 . The shooting requires a second boundary at η = η 0 , which we determine in the regime χ 1 as

Exact solution for α = π
In the special case that the wall opening angle is α = π , an exact solution for the flow depth can be obtained. We introduce the following conformal mapping of the physical space s + it = 2(x + iy) 1/2 (5.1) and in this case the domain now corresponds to s > 0. The governing equation (2.3) for the flow depth is recast as  This boundary condition suggests a solution of the form h = H(t) and upon substituting into (5.2), we find that H(t) satisfies the same ODE (4.2) and boundary conditions as χ (η). The complete solution is therefore given by where the argument of χ has positive sign when y > 0 and negative when y < 0. Contours of the height field (5.4) correspond to where χ (η * ) = h * . In particular the contact line y c (x) is given by η * = η 0 (see §4). We plot the contours (5.5) for h * = 0.01, 0.4, 0.7 in figure 4, showing excellent agreement with the numerically computed solution. We also note that when x 1,

The translated similarity solution
We can improve the agreement between the earlier similarity solution and the numerical computations for all α by translating it a distance x = −x 1 < 0 upstream. We use the similarity variable η = y/(x + x 1 ) 1/2 [19]. The idea is that this virtual origin accounts for the adjustment to the similarity solution in the region in which x ∼ 1. The solution χ (η) is unchanged but is translated upstream a distance x 1 . The parameter x 1 is chosen to agree with the contact line from the exact solution for α = π (5.4); we use x 1 = η 2 0 /4 ≈ 0.623. We first note that the translated similarity solution does not agree precisely with the exact solution for α = π away from the contact line because the contours take the form given by (5.5). Of course, the solution converges to the similarity shape for large x, as before.
The prediction for the contact line from the translated similarity solution is significantly improved from the untranslated similarity solution even for α = π (figure 5). The agreement for the contact line between the similarity solution and the numerical results could be improved further still by allowing the location of the virtual origin to depend on α. However, the single choice of x 1 = η 2 0 /4 provides sufficient agreement and enables the contact line to be simply approximated for any α.  (table 1). The prediction from the similarity solution (equation (6.2)) is plotted as a blue dashed line and we also include its behaviour as α → 0, which is d = (η 0 /α) 2

(a yellow dot-dashed line). (Online version in colour.)
Our prediction for the location of the contact line can be used to estimate the location of the detachment from the wall as a function of α. The similarity solution with virtual origin at −x 1 = −η 2 0 /4 predicts that the shape of the current edge where h = 0 is given by The intersection of the contact line (6.1) obtained from the similarity solution with the wall, y = x tan α, yields the following location for the detachment distance along the wall, This prediction is compared with the numerical results in figure 6. We note that locating the virtual origin at −x 1 = −η 2 0 /4 provides good agreement for all angles of expansion and the agreement is exact at α = π . As α → 0, equation (6.2) predicts that d ∼ (η 0 /α) 2 and we include this in figure 6. It provides good agreement over the whole range of opening angles, α.

Laboratory experiments
We conducted a series of experiments in which golden syrup was released from a line source at constant flux onto a Perspex tank inclined at an angle β = 10 • . The line source was restricted to one side of the tank, which was bounded by a wall parallel to the direction of steepest descent down the slope. Further downstream, the wall opened at its vertex with angle α, which was varied. The region around the wall was cleaned each time the angle was changed to ensure that the influence of surface tension at the detachment point was consistent across all experiments. The syrup was released from a lock gate behind which the depth was maintained fixed by hand to ensure constant flux. When the steady state was attained, the flow depth upstream, H ∞ , was obtained by projecting a narrow laser line onto the free-surface and recording the difference between the deflection of the line in the absence of fluid [12,20]. In order for this method to be effective the fluid must be relatively opaque and this required the addition of a few drops of paint to the syrup. It was confirmed (through rheometer measurements) that the syrup remained Newtonian with dynamic viscosity 89.6 Pa s.  The dimensional detachment distance along the wall, D, was measured from photographic images of the steady flow and the results are shown in table 1. They show good agreement with the numerical predictions and analytical estimates from our theory (figure 6).

Discussion and conclusion
This contribution has described the spreading of viscous fluid into an expanding channel using the equations of Stokes flow and the lubrication approximation. The governing equation has been solved numerically and far downstream the shape of the flow depth is self-similar, being governed by a balance between the gravity-driven flux down the slope and the cross-slope gradients in the hydrostatic pressure. We have obtained an exact analytical solution to the complete system in the case that there is a single, semi-infinite wall (α = π ). This solution motivates translating the similarity solution, which provides very good predictions for the location of the contact line and the wall detachment distance as a function of the vertex angle. These results also show good agreement with our laboratory experiments.
Our investigation was motivated by the need to optimize the design of lava flow barriers, which are used to protect towns and infrastructure in many regions of volcanic hazard. The flow around a vertex studied here is analogous to the flow at the downstream side of a deflecting barrier, albeit that our flows are uniform upstream. Civil defence authorities are particularly interested in the size and shape of the downstream dry region and our results demonstrate that surprisingly this may be insensitive to the shape of the trailing side of the barrier (the contact line is always well-approximated by equation (6.1)). This suggests that the size of the dry zone is primarily controlled by the cross-stream width of the barrier and the upstream deflection of the oncoming flow. Therefore, efforts should be directed at ensuring the barrier is as wide as possible in the cross-stream direction and shaped to minimize the upstream flow depth and potential for overtopping, a problem that was analysed in [12]. It would be interesting to extend the present analysis to downstream boundaries with curved shape and explore whether the current thickness is also insensitive to the shape in this case.