A solution method for one-dimensional shallow water equations using flux limiter centered scheme for open Venturi channels

Abstract

The one-dimensional shallow water equations were modified for a Venturi contraction and expansion in a rectangular open channel to achieve more accurate results than with the conventional one-dimensional shallow water equations. The wall-reflection pressure–force coming from the contraction and the expansion walls was added as a new term into the conventional shallow water equations. In the contraction region, the wall-reflection pressure–force acts opposite to the flow direction; in the expansion region, it acts with the flow direction. The total variation diminishing scheme and the explicit Runge–Kutta fourth-order method were used for solving the modified shallow water equations. The wall-reflection pressure–force effect was counted in the pure advection term, and it was considered for the calculations in each discretized cell face. The conventional shallow water equations produced an artificial flux due to the bottom width variation in the contraction and expansion regions. The modified shallow water equations can be used for both prismatic and nonprismatic channels. When applied to a prismatic channel, the equations become the conventional shallow water equations. The other advantage of the modified shallow water equations is their simplicity. The simulated results were validated with experimental results and three-dimensional computational fluid dynamics result. The modified shallow water equations well matched the experimental results in both unsteady and steady state.

Keywords

Shallow water equations wall-reflection pressure–force total variation diminishing scheme open Venturi channel contraction and expansion walls hydraulic jump flow depth

Introduction

The shallow water equations (SWEs) are used in various applications, such as river flow, dam break, open channel flow, etc. Compared to the 3D SWEs, 1D SWEs have a much lower cost in time-dependent simulations.¹ Kurganov et al.² introduced a semidiscrete central-upwind numerical scheme for solving the Saint-Venant equations, which is suitable for use with discontinuous bottom topographies.³ This scheme avoids the breakdown of numerical computation when the channel is at dry or near dry states. Another computational difficulty is that small flow depth leads to enormous velocity values near the dry states.³ By accurately calculating the wall-reflection pressure–force it is possible to prevent artificial acceleration of the flow.⁴ Spurious numerical waves propagate when the time discretization step is too large.^5,6 The total variation diminishing (TVD) method does not allow to increase total variation in time.⁷ According to Toro,⁷ the centered TVD scheme consists of a flux limiter blending of the FORCE scheme and the Richtmyer scheme. High-resolution schemes and flux limiters are suitable for avoiding phase error for monotone solutions.⁸ Partial differential equations can be solved by splitting them into a hyperbolic problem and a source problem.⁷ In the operator-splitting approach, the eigenvector projection and improved approaches are used for source term treatments.⁹ An open Venturi channel on the horizontal plane gives a subcritical flow regime before the Venturi contraction walls, and the flow regime changes from critical to supercritical after the Venturi expansion walls.¹⁰ As far as we know, a possibility of use of the flux limiter centered (FLIC) scheme for solving the 1D SWE at unsteady and steady states for nonprismatic channels is not available in the literature. The paper addresses this area with some modifications of SWE.

The underlying assumptions of 1D SWE are summarized as follows: velocity is uniform in the cross-section, water level in the cross-section is presented as a horizontal line, vertical acceleration is negligible, and streamline curvature is small. Therefore, pressure can be considered to be hydrostatic pressure.¹¹ Based on the conservation of mass and momentum, the Saint-Venant system of the conventional SWE can be written as^12,13

\frac{\partial A}{\partial t} = - \frac{\partial (A u)}{\partial x}

(1)

\frac{\partial (A u)}{\partial t} = - \frac{\partial (A u^{2})}{\partial x} - \frac{\partial (A h)}{\partial x} g + A g \sin α - A g S_{f}

(2)

Here, $A$ is the cross-sectional area perpendicular to the $x$ -direction, $u$ is the velocity, $h$ is the flow depth, $g$ is the acceleration due to gravity, $α$ is the channel inclination angle, and $S_{f}$ is the friction slope.

Here, modified SWEs are developed for the Venturi contraction and expansion for a rectangular channel. The centered TVD scheme is used for solving the modified SWE. MATLAB R2017a was used for the 1D simulations. The experiments were carried out in a trapezoidal open Venturi channel. The developed model is validated through experimental results without using analytical results. The paper proceeds from conventional SWE (“Modified 1-D SWEs for open Venturi channel” section) to the model development of modified SWE (“Centered TVD method for the modified 1D SWEs”) and the implementation of the TVD scheme for the modified SWE. The modified equations are then compared to the conventional SWE and validated with experimental results.

Modified 1D SWEs for open Venturi channel

A rectangular open Venturi channel is used for model development. The principle sketch is shown in Figure 1. The channel has a continuous bottom topography. In the Venturi section, the bottom width of the channel varies in the $x$ -direction. The walls are perpendicular to the bottom surface. The inlet of the channel is defined at $x = 0$ m. The walls are stationary. A control volume is shown in the Venturi contraction region.

Figure 1.

Principle sketch of the open Venturi channel with a rectangular cross-section. The selected control volume is in the Venturi contraction region. The top surface is open to the atmosphere.

Model development

Figure 2 shows a spatial and time discretized grid for one time step, which is based on the finite volume method.¹⁴ $n$ is the time index, $n \in {1, 2, \dots, N}$ . At time $t = n$ , $x$ coordinates are discretized. $j$ is the node index in the spatial grid, $j \in {1, 2, \dots, l}$ . We assume that $Δ x$ is a constant for all cells. The aim of an iteration step is to find the conserved variables at time $t = n + 1$ and $x = j$ from $t = n$ variable values. Time step is variable and is defined as $Δ t = t^{n + 1} - t^{n}$ .

Figure 2.

Semidiscretized grid, spatial discretization presented with $j$ -notations, and time discretization presented with $n$ notations.

The model emphasizes the impact of contraction and expansion walls. The control volume is in the Venturi contraction region to consider a maximum number of boundaries. Figure 3 shows the forces acting on the control volume in the $x$ -direction. The gravitational and the bottom wall-reflection forces act in the $z$ -direction and are not shown. If there is a channel downward inclination, a portion of the gravitational support is added in the $x$ -direction. The height of the control volume is equal to the flow depth, which varies from the inlet to the outlet of the control volume. Therefore, the inlet area ( $A_{j}$ ) and the outlet area ( $A_{j + 1}$ ) are not equal. Changes in flow depth from the inlet to the outlet are assumed to be linear at any given time. $R$ is the reflection pressure–force coming from the side boundaries. This component is not included in the conventional 1D SWE (see equation (2)). Sanders and Iahr⁴ noticed that the hydrostatic force coming from the channel walls in nonprismatic channels needs to be treated with precision, in order to avoid an artificial acceleration of the flow in the calculation. $R$ acts in the $y$ -direction only if there are no contraction or expansion regions. Therefore, the 1D conventional SWE can be applied to prismatic channels without any issue. Here, we assumed that the 1D conventional SWE does not consider the pressure–force coming from the sidewalls in $\frac{\partial (A h)}{\partial x} g$ . Easy numerical calculation and accurate pressure–force calculation are the advantages of having a separate term for the reflection pressure–force coming from the sidewalls. The wall-reflection pressure–force is equal to the hydrostatic pressure acting on the wall. $R$ is acting perpendicular to the sidewalls. According to the assumption that the changes in flow depth from the inlet to the outlet of the control volume are linear, we can assume that the wall-reflection pressure–forces coming from each of the two sidewalls of the control volume are equal. Here, $F_{f}$ is the resultant friction force from both sidewalls and the bottom surface. $β$ is the contraction angle of the control volume

\tan β = \frac{Δ b}{2 Δ x}

Δ b = b_{i - \frac{1}{2}} - b_{i + \frac{1}{2}}

Figure 3.

Forces on the control volume in the $x$ -direction, gravitational force is not shown.

$b$ is the bottom width of the channel. For a contraction and an expansion region, the sign of the $Δ b$ value becomes negative and positive, respectively. The sign of the $Δ b$ determines the sign of $R$ . $R$ acts opposite to the flow direction in a contraction region and in support of the flow direction in an expansion region.

Fundamental conservation laws are used for model development; the temperature is assumed to be constant at room temperature, and density is also assumed to be constant. Two conservation equations are produced by applying mass and momentum balances to the control volume. There is no difference in the mass balance equation compared to the conventional SWE, which is equal to equation (1). Applying the momentum balance to the control volume

\frac{\partial {(m u)}_{j}}{\partial t} = {(\dot{m} u)}_{j - \frac{1}{2}} - {(\dot{m} u)}_{j + \frac{1}{2}} + \sum_{i} F_{i}

(3)

Here, $m$ is the mass, $u$ is the velocity, and $F$ is the force. The $x$ -directional forces acting on the control volume are pressure–force $F_{p}$ , gravitational force $F_{g}$ , friction force $F_{f}$ , sidewall-reflection pressure–force $F_{R},$ and other external forces $F_{e}$ . Assuming there are no other external forces working on the control volume, $F_{e} = 0$ .

The pressure applied to the cross-section is a sum of atmospheric pressure and hydrostatic pressure. Atmospheric pressure is balanced from both sides of the cell faces as well as the top and bottom surfaces. The resultant hydrostatic pressure–force coming from the adjacent control volumes is

F_{p} = k_{g} ({(A h)}_{j - \frac{1}{2}} - {(A h)}_{j + \frac{1}{2}}) ρ g

$k_{g}$ is the ratio between the gravity height of the cross-sectional area to the flow depth. In this study, it is assumed that the gravity of the flow depth is half of its total height for a rectangular cross-sectional area. $F_{p}$ is not included as all pressure–forces apply on the control volume in the $x$ -direction. If the channel is a prismatic channel, $F_{p}$ includes all pressure–forces that act on the control volume in the $x$ -direction. Therefore, as explained above, the wall-reflection pressure–force coming from the sidewalls needs to be considered for nonprismatic channels

F_{R} = - 2 R \sin β = - 2 \sin β \int_{0}^{A_{s w, j}} p_{s w} d A

(4)

Here, $A_{s w}$ is the sidewall area of the control volume (one side). $p_{s w}$ is equal to the hydrostatic pressure acting on the sidewall. According to the assumptions, flow depth variation is linear from the inlet to the outlet of the control volume; therefore, average values can be used for the wall-reflection pressure–force and the area of the sidewalls, respectively

F_{R} \approx - 2 {\bar{p}}_{s w, j} {\bar{A}}_{s w, j} \sin β

(5)

The central differencing approach for the flow depth of a channel with a rectangular cross-section, $h_{j} = \frac{h_{j + \frac{1}{2}} + h_{j - \frac{1}{2}}}{2}$ , leads to

{\bar{A}}_{s w, j} = h_{j} \frac{Δ x}{\cos β}

{\bar{p}}_{s w, j} = k_{g} h_{j} ρ g

$F_{R}$ becomes

F_{R} = - 2 k_{g} Δ x h_{j}^{2} ρ g \tan β

(6)

When the channel has an inclined plane, a gravitational force acts with or against the flow direction. $α$ is the channel inclination angle. When the channel inclines downward, $α$ has a negative sign and the gravitational force acts with the direction of the flow. When the channel inclines upward, $α$ has a positive sign and the gravitational force acts against the direction of the flow. The sign of $F_{g}$ is decided by $α$ . $V$ is the volume of the control volume and is a function of $b$ , $h$ , and $x$ at any given time. The accurate volume calculation of the control volume is a very important step in identifying the different flow regimes in open nonprismatic channels.¹⁰ The friction slope $S_{f}$ gives the boundary friction force per unit weight of liquid present in the open channel.¹² According to Manning’s formula, the friction slope of the open channel can be presented as¹²

S_{f} = \frac{k_{M}^{2}}{k_{n}^{2} r_{h, j}^{4 / 3}} u_{j}^{2}

Here, $k_{M}$ is the Manning roughness factor, $r_{h}$ is the hydraulic radius, and $k_{n} = 1.0$ m^1/3 s⁻¹ is the unit corrector. For a rectangular channel, the hydraulic radius is $\frac{b h}{b + 2 h}$ . By substituting momentum and force terms into equation (3), the modified momentum balance equation can be stated as

\begin{matrix} \frac{\partial (A u)}{\partial t} = - \frac{\partial (A u^{2})}{\partial x} - k_{g} \frac{\partial (A h)}{\partial x} g \\ + k_{g} h^{2} g \frac{\partial b}{\partial x} + A g (\sin α - S_{f}) \end{matrix}

(7)

Compared to the conventional shallow water momentum balance equation, equation (2), the expression $2 k_{g} h^{2} g \frac{\partial b}{\partial x}$ is added to the new equation. The new term is related to the wall-reflection pressure–force in that it becomes zero when there is no contraction or expansion in the channel (when $\frac{\partial b}{\partial x} = 0$ ). In general, this term is only active in contraction or expansion regions.

Free falling at the end of the channel

In the experimental setup, the channel end was open, and the water was unhindered in flowing out of the channel. Accordingly, in the simulation, the physics of the last cells at the channel end needed to be modified with free falling properties. There is no friction effect when water does not touch the walls. Therefore, $S_{f} = 0$ , and a high gravitational force was added to the flow direction at the end of the channel. This effect is described by angle $α$ , and it reaches to $70 °$ – $90 °$ .

Centered TVD method for the modified 1D SWEs

The conservation equations are based on a rectangular open Venturi channel. The bottom width $b (x)$ is only a function of $x$ . The velocity, $u (x, b, t)$ and the flow depth $h (x, b, t)$ are functions of $x$ , $b (x),$ and time $t$ . The modified SWE (equations (1) and (3)) can be stated for a rectangular channel

\frac{\partial h}{\partial t} = - \frac{1}{b} \frac{\partial (h u b)}{\partial x}

(8)

\begin{array}{l} \frac{\partial (h u)}{\partial t} = - \frac{1}{b} \frac{\partial (b \frac{{(h u)}^{2}}{h})}{\partial x} - \frac{1}{2 b} \frac{\partial (b h^{2})}{\partial x} g \\ + \frac{1}{2 b} h^{2} g \frac{\partial b}{\partial x} + h g (\sin α - \frac{k_{M}^{2}}{k_{n}^{2}} {(\frac{b + 2 h}{b h})}^{\frac{4}{3}} {(\frac{h u}{h})}^{2}) \end{array}

(9)

In equation (9), $\frac{\partial b}{\partial x}$ is approximately equal to $2 \tan β$ , when $Δ x$ is not very small. It reduces the complexity of hyperbolic equation solving and allows for easy comparison with the conventional SWE. The compact form of the transport equations is

U_{t} + \frac{1}{b} F {(U)}_{x} = \frac{1}{b Δ x} S_{1} (U) + S (U)

(10)

Here

\begin{array}{l} U = (\begin{matrix} h \\ h u \end{matrix}) = (\begin{matrix} u_{1} \\ u_{2} \end{matrix}), \\ F (U) = (\begin{matrix} h u b \\ b \frac{{(h u)}^{2}}{h} + \frac{g b h^{2}}{2} \end{matrix}) = (\begin{matrix} b u_{2} \\ \frac{b u_{2}^{2}}{u_{1}} + \frac{g b u_{1}^{2}}{2} \end{matrix}), \\ S_{1} (U) = (\begin{matrix} 0 \\ h^{2} g \frac{Δ b}{2} \end{matrix}) = (\begin{matrix} 0 \\ \frac{Δ b}{2} g u_{1}^{2} \end{matrix}), \end{array}

\begin{matrix} S (U) = (\begin{matrix} 0 \\ h g (\sin α - \frac{k_{M}^{2}}{k_{n}^{2}} {(\frac{b + 2 h}{b h})}^{4 / 3} {(\frac{h u}{h})}^{2}) \end{matrix}) \\ = (\begin{matrix} 0 \\ u_{1} g (\sin α - \frac{k_{M}^{2}}{k_{n}^{2}} {(\frac{b + 2 u_{1}}{b u_{1}})}^{4 / 3} {(\frac{u_{2}}{u_{1}})}^{2}) \end{matrix}) \end{matrix}

(11)

Here, $U$ is the column vector of the conserved variables: $u_{1}$ and $u_{2}$ . $F (U)$ is the $x$ -directional column vector of the fluxes. $S (U)$ is the column vector of the source terms. The new term, which is the wall-reflection pressure–force, is presented as an extra source term, $S_{1}$ . The sidewall-reflection pressure–force effect is part of the pure advection term. However, it is written as a source term for the sake of easy numerical calculation. The pure advection term (advection flux and wall-reflection effect) is solved first with the centered TVD method. The source term is then solved with an ordinary differential equations (ODE) solver, the explicit Runge–Kutta fourth-order method. The initial condition for the ODE solver is the solution coming from the centered TVD scheme. Here, we emphasize the solving of the TVD scheme (not the ODE solving), because it is included in the wall-reflection pressure–force term. The variation of the bottom width ( $b$ ) presents the channel contraction and expansion effects numerically. The bottom-width variation effect is highlighted in the TVD solving method used here, compared to the conventional TVD solving method. The pure advection term is solved as

U_{j}^{n + 1} = U_{j}^{n} - \frac{Δ t}{b_{j} Δ x} [F {(U)}_{j + \frac{1}{2}}^{n} - F {(U)}_{j - \frac{1}{2}}^{n} - S_{1} {(U)}_{j}^{n}]

(12)

The source term is solved as

U_{j}^{n + 1} = U_{j, T V D}^{n + 1} + Δ t S (U_{j, T V D}^{n + 1})

(13)

Here, $U_{j, T V D}^{n + 1}$ is the solution coming from the TVD method. The FLIC scheme is used to calculate the flux of the cell face⁷

F {(U)}_{j + \frac{1}{2}}^{n} = F {(U)}_{j + \frac{1}{2}}^{n, L O} + \emptyset_{j + \frac{1}{2}} [F {(U)}_{j + \frac{1}{2}}^{n, H O} - F {(U)}_{j + \frac{1}{2}}^{n, L O}]

(14)

Here, $H O$ is the high-order flux, and $L O$ is the low-order monotone flux. The Richtmyer scheme ( $R I$ ) is used to calculate the high-order flux, and the FORCE scheme is used to calculate the low-order flux. $\emptyset$ is the flux limiter function. The high-order flux can be calculated as⁸

F {(U)}_{j + \frac{1}{2}}^{n, H O} = F (U_{j + \frac{1}{2}}^{n, R I})

(15)

U_{j + \frac{1}{2}}^{n, R I} = \frac{1}{2} [U_{j}^{n} + U_{j + 1}^{n}] + \frac{Δ t}{2 b_{j + \frac{1}{2}} Δ x} [F (U_{j}^{n}) - F (U_{j + 1}^{n})]

(16)

The low-order flux is based on the FORCE scheme⁷

F {(U)}_{j + \frac{1}{2}}^{n, L O} = F {(U)}_{j + \frac{1}{2}}^{n, F O R C E} = \frac{1}{2} [F {(U)}_{j + \frac{1}{2}}^{n, L F} + F (U_{j + \frac{1}{2}}^{n, R I})]

(17)

The Lax–Friedrichs ( $L F$ ) scheme can be presented as⁷

F {(U)}_{j + \frac{1}{2}}^{n, L F} = \frac{1}{2} [F (U_{j}^{n}) + F (U_{j + 1}^{n})] + \frac{b_{j + \frac{1}{2}} Δ x}{2 Δ t} [U_{j}^{n} - U_{j + 1}^{n}]

(18)

The SUPER-BEE flux limiter ( $\emptyset$ ) is used as the flux limiter.^7,15 The cell face flux depends on the value of flow parameter $r$ , which is the flux limiter function. The general idea of $r$ is that it is the ratio of upwind change to local change. Here, a total energy-based method is used to calculate these changes. $q$ is defined as the total energy per unit mass which is the sum of potential and kinetic energies per unit mass

q = \frac{1}{2} h g + \frac{1}{2} u^{2} = \frac{1}{2} u_{1} g + \frac{1}{2} {(\frac{u_{2}}{u_{1}})}^{2}

(19)

Time step is related to wave propagation speed. $S_{m a x}$ is the maximum wave propagation speed. The wave propagation speed is calculated from the Froude number. The time step can be calculated as

S^{n} {m a x}_{=} m a x (a b s (u_{j} + \sqrt{h_{j} g})), \forall j

(20)

Δ t = \frac{c_{m a x} Δ x}{S^{n}_{m a x}}

(21)

According to the observations, Courant numbers higher than 0.7 led to high numerical diffusions at $Δ x = 0.01$ m. A constant time step gave very diffusive and inaccurate results; therefore, a variable time step was used instead.

Modified versus conventional SWEs for open Venturi channels

To compare the advantages of the modified SWE over conventional SWE, measurements and data from the experimental setup are used to supply the necessary variables. The results validate the modification. In the next section, the calculated results will further be compared to the measured results of the experiment as well as to the modeled results of 3D computational fluid dynamics (CFD). A more detailed description of the experimental setup will be given in section 6; here we only consider the values necessary for the calculations.

The total length of the channel is 3.7 m. The Venturi contraction region is $2.95 m < x < 3.1 m$ , and the expansion region is $3.3 m < x < 3.45 m$ . The width of the channel along the $x$ -axis is given as

0 m \leq x \leq 2.95 m : b = 0.2 m,

2.95 m < x < 3.1 m : b = 0.2 - \frac{x - 2.95}{1.5},

3.1 m \leq x \leq 3.3 m : b = 0.1 m,

3.3 m < x < 3.45 m : b = 0.1 + \frac{x - 3.3}{1.5},

3.45 m \leq x \leq 3.7 m : b = 0.2 m

At initial conditions, the whole channel was filled with water, and all the node points were measuring the same flow depth and zero velocity. According to this condition, there was no flux propagation. However, in the contraction and expansion regions, the conventional SWE produced a flux difference, because of $\frac{\partial (A h)}{\partial x} g = h^{2} g \frac{\partial b}{\partial x}$ for constant $h$ . At constant flow depth and zero velocity, the expression $\frac{\partial b}{\partial x}$ can erroneously produce a flux gradient in a contraction or an expansion region. This is because the conventional SWE do not account for the wall-reflection pressure–force effect. In the modified SWE, this error is avoided by considering the wall-reflection pressure–force effect in the pure advection term.

Figure 4 shows a comparison between conventional SWE and modified SWE results after 0.01 s. The initial flow depth was 100 mm in the whole channel with zero velocity. At the contraction and expansion regions, flow depths and velocities change considerably in the conventional SWE. Moreover, these variations expand with time and extend into the whole channel. This error produces inaccurate results.

Figure 4.

Comparison between conventional and modified SWE with zero velocity and constant flow depth in the whole channel at initial conditions. Results after 0.01 s: (a) flow depth along the $x$ - axis and (b) velocity along the $x$ - axis.

Model validation

The calculated results for the modified SWE are further validated by experimental and 3D CFD results.

Model validation with experimental results

The open Venturi channel used for the experimental model validation was located at University College of Southeast Norway. Level transmitters were located along the central axis of the channel. The transmitters were movable along the axis. The accuracy of the Rosemount ultrasonic 3107 level transmitters was ±2.5 mm for a measured distance of less than 1 m.¹⁰ The channel had a trapezoidal shape with a trapezoidal angle of $70 °$ . In modeling and simulations, a rectangular channel was used with all the other dimensions equal to the experimental setup. All of the experimental values presented in this paper are average values of sensor readings taken over a period of 5 min at each measuring point. It is possible to change the channel inclination angle ( $α$ ); negative values for $α$ indicate a downward inclination.

From the inlet, water was added to the channel at the horizontal plane ( $α = 0 °$ ) until the flow reached quasi-steady state while the outlet was kept open. Figures 5 and 6 show the comparison of simulated and experimental results. The constant flow rate at the inlet was 6.67 kg/s. Flow depth and velocity at the inlet were constant values of, respectively, 100 mm and 0.33 m/s, which were approximately equal to the inlet conditions of the experimental setup. The last five computation cells had free falling properties to account for the free falling water at the end of the channel of the experimental setup. The wall-reflection pressure–force coming from the Venturi contraction walls could be clearly observed while the channel was at the horizontal plane. In this condition, the gravitational force is zero. When water collapsed onto the Venturi contraction walls, a hydraulic jump was propagated. The hydraulic jump travelled upstream until it reached the inlet. With the channel at the horizontal plane, this hydraulic jump propagation was caused only by the wall-reflection pressure–force and friction–force effects.^10,16 According to the critical depth calculation in Welahettige et al.,¹⁰ the average critical flow depth before the contraction region is 48.4 mm. Figure 5 shows a dynamic situation of the flow profile. The contraction wall-reflection pressure–force acts opposite to the flow direction. This is the reason for the hydraulic jump travelling upstream. The wall-reflection pressure–force causes the flow regime to change from supercritical flow to subcritical flow in the hydraulic jump.¹⁶ The wall-reflection pressure–force coming from the contraction walls reduces the speed of the upstream flow, which causes the change in the flow regime.

Figure 5.

The hydraulic jump is moving upstream due to the wall-reflection pressure–force effect coming from the Venturi contraction walls, a dynamic result after 8.9 s. The arrows show the traveling direction of the hydraulic jump. The flow direction is opposite to the direction of the arrow. (a) The simulated flow depth result along the central axis and (b) the experimental results.

Figure 6.

Quasi-steady state results, flow rate at 400 kg/min; (a) The simulated flow depth results along the central axis, (b) the simulated velocity results along the central axis, and (c) the experimental result.

Figure 6 shows the quasi-steady state results, following Figure 5. The flow depth comparison between the simulated and the experimental results in Figure 6(a) indicates the accuracy of the modified SWE. The modified SWE result is well matched with the experimental results compared to the conventional SWE. The channel at the horizontal plane, the Venturi contraction can cause a significant change in the flow regime. According to the flux calculation in the contraction and the expansion regions, the results from the conventional SWE deviated from the experimental results at quasi-steady state. The wall-reflection pressure–force coming from the contraction walls changed the flow regime from supercritical to subcritical, whole channel section before the contraction region. The velocity profile in Figure 6(b) can be used to explain the wall-reflection pressure–force effect in the Venturi expansion region. At the expansion region ( $3.3 < x < 3.45$ m), the velocity drastically increased due to the wall-reflection pressure–force effect. At the channel expansion, the wall-reflection pressure–force effect comes in support of the flow direction, which results in increased velocity. Before the contraction region ( $x < 2.95$ m), the velocity reduced compared to the expansion region. This indicates the sign of the wall-reflection pressure–force: at the Venturi contraction region the sign is negative, and at the Venturi expansion region it is positive.

Model validation with 3D CFD result

Further, the result of the modified SWE result was compared to a CFD result.¹⁶ Three-dimensional CFD simulations are based on the volume of fluid method. Water and air are the materials in the fluid domain. An artificial compression term is activated at the interface.^17–19 Time discretization is based on the implicit Euler method for a transient calculation. Pressure–velocity coupling is based on the SIMPLE scheme with a second-order upwind correction. The standard k–ε model is used for the turbulence handling. Wall surface roughness is used to calculate the wall friction, which is 15 µm. The mesh contains 0.74 million elements with a maximum cell size of 10 mm. ANSYS Fluent 16.2 (commercial code) was used as the simulation tool.^10,16 The 3D CFD study was done with the same experimental setup for 100 kg/min flow rate. The quasi-steady state results are shown in Figure 7. The modified SWE result was well matched with the CFD result compared to the conventional SWE. A similar flow profile was achieved by Berg et al.²⁰ from CFD simulation for an open Venturi channel. The error propagated in the contraction and expansion regions caused the high deviation when using the conventional SWE.

Figure 7.

A comparison with 3D CFD result with the modified and the conventional shallow water results, quasi-steady state, the flow rate at 100 kg/min. The CFD result is from Welahettige et al.¹⁶ CFD: computational fluid dynamics.

Figure 8 shows a flow surface for the full channel (iso-surface of water volume fraction of 0.5) from the 3D CFD simulations,¹⁰ which is related to Figure 6(a). The flow rate is 400 kg/min, and the channel is at a horizontal angle. One-dimensional simulation surface profile is well matched with the 3D CFD surface profile.

Figure 8.

Water flow rate 400 kg/min and open channel at horizontal position: Simulated flow surface for the full channel (iso-surface of water volume fraction of 0.5). The flow direction is left to right.¹⁰

Conclusion and future work

The 1D conventional SWEs cannot be applied to a channel with a contraction and an expansion region (Venturi channel). Because conventional SWE neglect the wall-reflection pressure–force effect, they are suitable for prismatic channels only. The modified 1D SWEs are developed by considering the wall-reflection pressure–force effect. The modified SWEs can be applied to both prismatic and nonprismatic channels, especially those with contraction and the expansion regions.

This study will further extend into drilling fluid flow measurement in an open Venturi channel. The non-Newtonian properties of drilling fluid will be considered. Further, the scenario of a reflection hydraulic jump hitting the inlet will be considered in the future study.

Footnotes

Acknowledgment

The authors gratefully acknowledge the resources for experiments and simulations provided by the University College of Southeast Norway.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Economic support from The Research Council of Norway and Statoil ASA through project no. 255348/E30 “Sensors and models for improved kick/loss detection in drilling (Semi-kidd)” is gratefully acknowledged.

References

Vreugdenhil

CB.

Numerical methods for shallow-water flow. Dordrecht: Springer Netherlands, 1994.

Kurganov

Noelle

Petrova

et al . Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton–Jacobi equations. SIAM J Sci Comput 2001; 23: 707–740.

Kurganov

Petrova

A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system. Commun Math Sci 2007; 5: 133–160.

Sanders

Iahr

High-resolution and non-oscillatory solution of the St. Venant equations in non-rectangular and non-prismatic channels. J Hydraul Res 2001; 39: 321–330.

Bermudez

Vazquez

ME.

Upwind methods for hyperbolic conservation laws with source terms. Comput Fluids 1994; 23: 1049–1071.

Vázquez-Cendón

ME.

Improved treatment of source terms in upwind schemes for the shallow water equations in channels with irregular geometry. J Comput Phys 1999; 148: 497–526.

Toro

EF.

Riemann solvers and numerical methods for fluid dynamics – a practical introduction. 3rd ed. Heidelberg: Springer Science & Business Media, 2009.

LeVeque

RJ.

Finite volume methods for hyperbolic problems. 1st ed. Cambridge: Cambridge University Press, 2002.

Tseng

M-H.

Improved treatment of source terms in TVD scheme for shallow water equations. Adv Water Resour 2004; 27: 617–624.

10.

Welahettige

Lie

Vaagsaether

Flow regime changes at hydraulic jumps in an open Venturi channel for Newtonian fluid. J Comput Multiphase Flows 2017; 9(4): 169–179.

11.

Aldrighetti

Computational hydraulic techniques for the Saint Venant equations in arbitrarily shaped geometry. University of Trento, 2007. http://eprints.biblio.unitn.it/1395/1/PhDTS_n.52_.pdf

12.

Akan

AO.

Open channel hydraulics. 1st ed. Burlington: Elsevier/BH, 2006.

13.

Szymkiewicz

Numerical modeling in open channel hydraulics. New York: Springer, 2010.

14.

Versteeg

Malalasekera

An introduction to computational fluid dynamics: the finite volume method. 2nd ed. Gosport, Hants: Pearson Education Ltd, 2007.

15.

Roe

PL.

Characteristic-based schemes for the Euler equations. Annu Rev Fluid Mech 1986; 18: 337–365.

16.

Welahettige

Lie

Vaagsaether

Computational fluid dynamics study of flow depth in an open Venturi channel for Newtonian fluid. In: Proceedings of the 58th SIMS, pp.29–34. Reykjavik: Linköping University Electronic Press.

17.

Ubbink

Numerical prediction of two fluid systems with sharp interfaces. PhD Thesis, University of London, UK, 1997.

18.

Rusche

Computational fluid dynamics of dispersed two-phase flows at high phase fractions. Exhibition Road, London: Imperial College London (University of London), 2003.

19.

Weller

Tabor

Jasak

et al . A tensorial approach to computational continuum mechanics using object-oriented techniques. Comput Phys 1998; 12: 620–631.

20.

Berg

Malagalage

Agu

et al . Model-based drilling fluid flow rate estimation using Venturi flume. IFAC-PapersOnLine 2015; 48: 171–176.