Next Article in Journal
Does Rural Water System Design Matter? A Study of Productive Use of Water in Rural Nepal
Next Article in Special Issue
Investigation on the Effect of Structural Parameters on Cavitation Characteristics for the Venturi Tube Using the CFD Method
Previous Article in Journal
Feasibility Assessment of a Water Supply Reliability Index for Water Resources Project Planning and Evaluation
Previous Article in Special Issue
Wall Stresses in Cylinder of Stationary Piped Carriage Using COMSOL Multiphysics
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Technical Note

A Simple and Unified Linear Solver for Free-Surface and Pressurized Mixed Flows in Hydraulic Systems

1
School of Hydropower and Information Engineering, Huazhong University of Science and Technology, Wuhan 430074, China
2
Yangtze River Scientific Research Institute, Wuhan 430010, China
*
Author to whom correspondence should be addressed.
Water 2019, 11(10), 1979; https://doi.org/10.3390/w11101979
Submission received: 19 August 2019 / Revised: 18 September 2019 / Accepted: 19 September 2019 / Published: 23 September 2019
(This article belongs to the Special Issue Pipeline Fluid Mechanics)

Abstract

:
A semi–implicit numerical model with a linear solver is proposed for the free-surface and pressurized mixed flows in hydraulic systems. It solves the two flow regimes within a unified formulation, and is much simpler than existing similar models for mixed flows. Using a local linearization and an Eulerian–Lagrangian method, the new model only needs to solve a tridiagonal linear system (arising from velocity-pressure coupling) and is free of iterations. The model is tested using various types of mixed flows, where the simulation results agree with analytical solutions, experiment data and the results reported by former researchers. Sensitivity studies of grid scales and time steps are both performed, where a common grid scale provides grid-independent results and a common time step provides time-step-independent results. Moreover, the model is revealed to achieve stable and accurate simulations at large time steps for which the CFL is greater than 1. In simulations of a challenging case (mixed flows characterized by frequent flow-regime conversions and a closed pipe with wide-top cross-sections), an artificial slot (A-slot) technique is proposed to cope with possible instabilities related to the discontinuous main-diagonal coefficients of the linear system. In this test, a slot-width sensitivity study is also performed, and the suitable slot-width ratio (ε) for the linear solver is suggested to be 0.05–0.1.

1. Introduction

Hydraulic systems are often filled with free-surface and pressurized mixed flows. For example, the urban drainage systems are generally designed to convey water under free-surface conditions. However, during intense rain events, water may overflow from storm sewers through manholes onto the land surface and then cause inundation. In this case, the sewer is often filled with mixed flows (characterized by the simultaneous occurrence of free-surface and pressurized flows). Mixed flows are highly dynamic [1], with frequent transitions of flow regimes. Moreover, characteristics of the two regimes are quite different, e.g., the celerity of acoustic waves of pressurized flows may be 2–3 orders of magnitude larger than that of free-surface flow waves. In practice, the fluid dynamics of the mixed flows are usually simulated using one-dimensional (1D) hydrodynamic models, where accurate and fast simulations are expected to provide real-time data for storm management. However, a unified mathematical description of the free-surface and the pressurized parts of the mixed flows, as well as a simple solution, are often challenging [2].
The existing models for mixed flows can be divided into two categories, which, respectively, describe the flows using two sets of equations (TSE) and a single set of equations (ASE). A TSE model solves the free-surface and the pressurized parts of mixed flows separately, mainly including the rigid column approach [3,4,5] and the shock-fitting approach [6,7,8,9,10,11]. This kind of model explicitly tracks the pressurization front, and their advantages and limitations have been well reviewed in Bousso et al. [12]. Using a single set of equations, an ASE model builds a unified formulation applied to different regimes of the mixed flow. In ASE models, the Saint–Venant equations (SVE) or their equivalents are often used as the mathematical descriptions. In terms of the form of the ASE, two sub kinds of these models (ASE-1 and ASE-2) have been identified in references.
In ASE-1 models, the pressure gradient in the momentum equation is divided into a gradient of water depth and a source term associated with bed slope, while the governing equations are written in their conservative form [2,13,14,15,16,17,18]. The ASE-1 models often adopt Godunov-type schemes and Riemann solvers (e.g., Bourdarias and Gerbi [19]), so they are good at capturing physical discontinuousness and are able to automatically track the pressurization fronts in mixed flows. As a result, the ASE-1 models are often called the “shock-capturing” family of models. However, the governing equations used in ASE-1 models may not be directly used to simulate pressurized flows, because the free surface term is not explicitly included. Hence, additional techniques are needed for the ASE-1 models to quantify the piezometric head when the pipe is pressurized. Two popular techniques are the Preissmann slot technique [20] and the two-component pressure approach [17,21]. The first technique assumes a narrow slot is in the ceiling of the pipe, while the second decouples the hydrostatic pressure from surcharged pressures.
In ASE-2 models, the pressure term is explicitly expressed in the momentum equation instead of the two terms of water depth and bed slope, and a nonconservative form of equations is used [22,23,24]. The pressure represents either the water level for free-surface parts or the piezometric head for pressurized parts, so a unified description of mixed flows is reached. Moreover, benefiting from the nonconservative form of momentum equations, the advection term can be calculated explicitly using an Eulerian–Lagrangian method (ELM). Using such kind of governing equations, Casulli and Stelling [23] proposed a semi–implicit formulation that naturally applies to the two regimes of mixed flows, where the pressure (water levels for free-surface parts and piezometric heads for pressurized parts) is synchronously solved by a coupling solution of the continuity and momentum equations. Because the cross-sectional conveyance area is a nonlinear function of water level, the velocity–pressure coupling in their model results in a nonlinear system which requires iterative solutions. Benefiting the high stability and large time steps, the ASE-2 model is reported to be tens of times faster than conventional models for mix flows in [25].
On the other hand, the linearization technique (based on the solutions of time level n) is often used in existing models for mixed flows, and should be noticed. For example, in Fuamba [9] and SWMM [26], the A / t (A is cross-sectional conveyance area) in the continuity equation is linearized. In Ji [15], all nonlinear coefficients for governing variables are linearized based on the solutions of time level n. In Bourdarias and Gerbi [19], the numerical flux at cell-faces and the source term are linearized using the solutions of time level n, to construct a linear system. The numerical formulation and solution are sharply simplified by various linearization techniques in these models.
In this study, using a linearization [9,26], a semi–implicit numerical model with a linear solver is proposed for the free-surface and pressurized mixed flows in hydraulic systems, which solves the two flow regimes within a unified formulation and very simple. The new model may be considered as a variant of the ASE-2 model of [23]. In the new model, the velocity–pressure coupling is solved by constructing and solving a linear system instead of the nonlinear system in [23]. Characteristics of representative types of mixed flows in pipes are analyzed, and the applicability of the new model to them is discussed. The model is then tested using three classical problems of mixed flows.

2. Methods

The governing equations, the description of mixed flows, and the numerical formulation of the proposed new model for mixed flows are introduced.

2.1. Governing Equations

The interest of this paper is confined to 1D modeling of incompressible mixed flows in closed inelastic pipes and connected open channels. Assuming hydrostatic flow, the governing equations (including the continuity and momentum equations) for such mixed flows are taken to be:
A t + A u x = 0
u t + u u x = g η x g n m 2 u | u | R 4 / 3
where u = averaged velocity of cross-section; A = wetted cross-sectional area (conveyance area); x = longitudinal distance along the channel; t = time; η = the pressure representing either the water level (for free-surface parts) or the piezometric head (for pressurized parts); g = gravity acceleration; nm = Manning’s roughness coefficient; and R = hydraulic radius.
According to Casulli and Stelling [23], using Equations (1) and (2), a unified description of the free-surface parts and the pressurized parts of mixed flows is reached. In the free-surface parts of mixed flows, the pressure η represents the water level, and the pressure varies synchronously with respect to the water level. In the pressurized parts of mixed flows, the pressure η represents the piezometric head which includes two parts, respectively, the hydrostatic part and the surcharged part. In fact, when the flow in a pipe is pressurized, the hydrostatic part can be determined using the information of height and geometry of the pipe, after the water compressibility and wall elasticity have been neglected. As a result, only the surcharged part of the piezometric head is left unknown. In [23] and this study, for the pressurized parts of the mixed flows, the piezometric head is used directly (namely the pressure η is not further divided) for simplicity.
Using the unified description of the mixed flows based on Equations (1) and (2), the pressure (including the water levels for free-surface parts and the piezometric heads for pressurized parts) is universally quantified by solving the velocity-pressure coupling, in the nonlinear solver of Casulli and Stelling [23]. In the new model, this idea of [23] is also adopted, but the velocity–pressure coupling is solved by constructing and solving a linear system instead of a nonlinear system.
The A in Equation (1) is a nonlinear function of water levels. If the continuity equation with a nonlinear A is directly coupled with the momentum equation, the resulting system is to be nonlinear and needs to be solved iteratively. To avoid the complex solutions of a nonlinear system, the term A / t in Equation (1) is linearized in this study in a similar way of Fuamba [9] and SWMM [26]. The resulting continuity equation after a local linearization is given by:
B η t + A u x = 0
where B (= A / η ) represents the wet surface width of a cross-section for the free-surface parts and is equal to 0 for the pressurized parts of a closed pipe. In Equations (2) and (3), for a cross-section located in a closed pipe, A is a nonnegative nonlinear function of the water level when the pipe is not fully filled; A is equal to the area of the cross-section (Ap) when the water level is equal to or larger than the celling of the cross-section (zc).
It must be pointed out that: the local linearization of the continuity equation does not break the unified description and solution of the mixed flows using Equations (1) and (2). Based on Equations (2) and (3), the pressures of the free-surface and the pressurized parts of the mixed flows can also be universally quantified by solving the velocity-pressure coupling.

2.2. Computational Grid and Description of Mixed Flows

A channel with both closed and nonclosed reaches is considered here, which is divided by a set of 1D cells (elements). The ne and ns are used to denote, respectively, the quantities of cells and sides in the computational grid, and every cell/side has a unique global cell/side index. A staggered grid of variable arrangement is used. The velocity u is defined at side centers (i − 1/2, i + 1/2, …) and the pressure η is defined at cell centers (i − 1, i + 1, …), as shown in Figure 1. A cross-section with general topographical data is arranged at the center of each cell, and the geometry of a cell is described by the following variables. For cell i, Δxi represents the length of the cell; Bi is the wet surface width and is related to the cell water level that varies with respect to time; Ai represents the conveyance area of the cross-section located at the center of the cell. For side i + 1/2 which connects cells i and i + 1, Δxi+1/2 is the distance between the centers of cells i and i + 1; Ai+1/2 and Ri+1/2 are, respectively, the conveyance area and the hydraulic radius of side i + 1/2, which are interpolated using those of cells i and i + 1.
The cells in closed and nonclosed reaches are distinguished by a mark (CLO). CLO is set to 1 for cells in closed reaches and set to 0 for those in nonclosed reaches. In closed reaches, cells are further distinguished according to their flow regimes (denoted by a mark “PRE”). PRE is set to 1 for cells in the pressurized parts of the closed pipe reach and set to 0 for those with a free-surface flow regime (see Figure 1). The cells in all the free-surface and the pressurized parts of the channel together contribute to the construction of an algebraic system.

2.3. Numerical Formulation

2.3.1. Discretization of Governing Equations

The momentum equation is discretized using an operator-splitting technique, and the model solves the governing equations in three steps. First, all the explicit terms in the momentum equation are computed to obtain provisional velocities. Second, the velocity–pressure coupling is solved to obtain new pressures. Third, a back substitution of the new pressures into the momentum equation is performed to obtain the final velocities.
The advection term in the momentum equation is solved by a pointwise ELM. The ELM is an explicit advection scheme, combining the stability and accuracy of Lagrangian approaches with the convenience of the fixed computation grid of Eulerian approaches. The inherently upwind property of the ELM eliminates the Courant–Friedrichs–Lewy number (CFL) restriction associated with explicit Eulerian methods and allows large time steps. Solution of the advection term using an ELM requires a backtracking of cell faces. Because the velocity varies along a channel (as a result of the varying cross-sectional topography, the friction, etc.), a multistep tracking technique [27] is used to ensure the accuracy of the ELM tracking.
The backtracking of the ELM starts at a side center, and a multistep backward Euler technique is performed to backtrack from the known location at time level n + 1 to the unknown foot of the backward trajectory at time level n. The tracking displacement of each substep is given by Δ x = u Δ t / N b t , where “Nbt” is the number of substeps and u is the velocity vector at the starting point of each substep. When a tracking point has been found, the velocity at this point is interpolated based on the velocity field of time level n, for calculating the next displacement and searching the next tracking point. Once the location of the foot of the trajectory at time level n is found, the initial condition at that time step is interpolated and recorded as ubt. The cell-face velocity is then directly updated with ubt, which means that the advection term is solved.
To eliminate the stability restrictions imposed by rapid gravity waves, the θ semi-implicit method [28,29,30,31] is used to discretize the gradient of the pressure in two parts, where the explicit and the implicit parts have a weight factor of 1 − θ and θ, respectively. For both the two regimes of mixed flows, the momentum equation can be universally discretized as (side i + 1/2):
( 1 + g Δ t n m 2 | u b t , i + 1 / 2 n | R i + 1 / 2 n 4 / 3 ) u i + 1 / 2 n + 1 = u b t , i + 1 / 2 n ( 1 θ ) g Δ t η i + 1 n η i n Δ x i + 1 / 2 θ g Δ t η i + 1 n + 1 η i n + 1 Δ x i + 1 / 2
where θ is the implicit factor; Δt is the time step; superscript n indicates time level n; and the riverbed friction term is discretized using ubt and un+1 to enhance computational stability.
When all the explicit terms in the discretized momentum equations are incorporated, the unknowns (the pressures at time level n + 1, ηn+1) emerge and Equation (4) is then transformed into the following form:
u i + 1 / 2 n + 1 = G i + 1 / 2 n F i + 1 / 2 n θ Δ t g 1 F i + 1 / 2 n η i + 1 n + 1 η i n + 1 Δ x i + 1 / 2
where F i + 1 / 2 n = 1 + g Δ t n m 2 | u b t , i + 1 / 2 n | / R i + 1 / 2 n 4 / 3 ; G i + 1 / 2 n is the incorporated explicit terms and G i + 1 / 2 n / F i + 1 / 2 n can be regarded as the provisional solution by solving all the explicit terms, where G i + 1 / 2 n = u b t , i + 1 2 n ( 1 θ ) g Δ t η i + 1 n η i n Δ x i + 1 / 2 .
Before the introduction of the discretization of the continuity equation, we want to explain that the different derivates (related to the term A in Equation (1)) have different roles in terms of their physics meanings. The term A / t describes the change of the cross-sectional area (or the water volume of a cell) with respect to time. Hence, in the discretization of the term A / t , the A i n + 1 and the A i n of a cell, respectively, at time levels n and n + 1 should be both involved. When the A / t is transformed into the B η / t in Equation (3), it means that η i n + 1 and η i n of a cell should be both involved in the discretization of the B η / t . At the same time, in the discrete continuity equation, the resulting discretized expression of the term ( A u ) / x accounts for the calculation of the water fluxes through cell faces. For the discretization of the ( A u ) / x , we can use either the A i ± 1 / 2 n or the A i ± 1 / 2 n + 1 (at cell faces), which may only influence the order of the accuracy of solutions. In practice, the A i ± 1 / 2 n is often used for simplicity (such as [23]), and it is also used in our model.
The continuity equation is discretized using a finite-volume method, where the mass of water is conserved both locally and globally by the unique flux at a cell face. For cell i, the discrete pressure ηi is assumed to be constant within the cell. Moreover, the wet surface width of cell i (Bi) is explicitly represented by the value of time level n ( B i n ) within a time step (Δt). From time level n to n + 1, a finite-volume discretization of the continuity equation, Equation (3), is given by:
B i n Δ x i ( η i n + 1 η i n ) = θ Δ t ( A i + 1 / 2 n u i + 1 / 2 n + 1 A i 1 / 2 n u i 1 / 2 n + 1 ) ( 1 θ ) Δ t ( A i + 1 / 2 n u i + 1 / 2 n A i 1 / 2 n u i 1 / 2 n )
The nonlinear terms (Bi and Ai±1/2) are present, respectively, on the left hand side (LHS) and the right hand side (RHS) of Equation (6). However, when the Bi on the LHS and the Ai±1/2 on the RHS have been explicitly represented by B i n and A i ± 1 / 2 n (within the Δt, from time level n to n + 1), the solution problem associated with the nonlinear terms does not exist any longer.
Our linear solver and the nonlinear solver in [23] use different discrete continuity equations. In [23], the LHS of the discrete continuity equation may written as Δ x i [ A i ( η i n + 1 ) A i ( η i n ) ] , where the function A(η) is a nonlinear function of water level. In the LHS of the discrete continuity equation in [23], existence of the nonlinear function A(η) potentially leads to a complex formulation for the construction and solution of the velocity-pressure coupling.

2.3.2. Solution of Velocity–Pressure Coupling

With the advection term explicitly calculated (by the ELM) and the A / t linearized, the coupling solution of the continuity and momentum equations in the new model only needs to solve a linear system. The velocity–pressure coupling is performed by substituting Equation (5) into Equation (6). The resulting equation is given by (i = 1, 2, …, ne):
C i a η i 1 n + 1 + C i b η i n + 1 + C i c η i + 1 n + 1 = B i n Δ x i η i n + r i n
where C i a = g θ 2 Δ t 2 A i 1 / 2 n F i 1 / 2 n Δ x i 1 / 2 , C i c = g θ 2 Δ t 2 A i + 1 / 2 n F i + 1 / 2 n Δ x i + 1 / 2 , C i b = B i n Δ x i C i a C i c , r i n = θ Δ t [ A i + 1 / 2 n G i + 1 / 2 n / F i + 1 / 2 n A i 1 / 2 n G i 1 / 2 n / F i 1 / 2 n ] ( 1 θ ) Δ t [ A i + 1 / 2 n u i + 1 / 2 n A i 1 / 2 n u i 1 / 2 n ] .
Equation (7) forms a tridiagonal linear system of ne equations with the cell pressures (ηn+1) as unknowns. When all the dry cells have been excluded from the system, for wet cells we have C i a < 0 , C i c < 0 and C i b > 0 . As a result, this linear system has a symmetric and positive definite coefficient matrix, and can be solved using a direct solution method without iterations, to obtain the final solutions (ηn+1). Once the ηn+1 have been obtained, the velocity un+1 defined at a side center is calculated by Equation (5), and then the Ai and the Ri are updated.
Relative to the solution of a linear system, the solution of a nonlinear system (e.g., [23]) which requires iterations will be much more complex. Casulli and Stelling [23] found that the convergence of most Newton-type methods is often problematic for the solution of mixed flows, and they used a nested Newton-type method to solve their nonlinear system.
Equation (7) is applicable to both free-surface and pressurized flow regimes, based on which a 1D unified numerical model is developed for mixed flows. However, in preliminary tests of certain mixed flows, it is found that the main-diagonal coefficients of the matrix of the linear solver may be discontinuous in value, arising from the potentially abrupt change of wet surface widths of cross-sections. The problem may leads to nonphysical oscillations in the solutions, or even unstable simulations. In the next section, characteristics of several typical mixed flows are analyzed, and the applicability of the new model to each type of them is discussed.

2.3.3. Considerations of Abrupt Change of Wet Surface Width in Mixed Flows

The wet surface width of a cross-section (Bi) has a positive finite value for free-surface flows, and is equal to 0 for pressurized flows in pipes. Hence, abrupt changes of the Bi may occur in mixed flows, especially when the cross-section of a closed pipe has a wide top (e.g., the rectangular cross-section in Figure 2). In Figure 2, the Bi changes abruptly from the width of Wi (representative width of the cross-section) in the free-surface part to 0 in the pressurized part, at the transition point. However, if the pipe has a closed circular cross-section, the Bi is a decreasing function of the water level and gradually approaches 0 when the water level approaches the ceiling of the pipe. Such a cross-section shape is defined as the “shrinking-top” here. In case of a shrinking-top cross-section, variation of the Bi shows a gradual change near the pipe ceiling instead of an abrupt change, at a transition point of different flow regimes (see Figure 1).
When the current model is used to simulate pure free-surface flows, the Bi has a positive value and the coefficient matrix of the linear system arising from the velocity–pressure coupling is diagonally dominant. The larger the Bi is, the larger the main-diagonal coefficient is. When the model is used to simulate pure pressurized flows, with Bi = 0, the coefficient matrix is no more diagonally dominant. When the model is used to simulate mixed flows in a closed pipe with wide-top cross-sections, abrupt changes of the Bi at a transition point will result in a linear system with discontinuous main-diagonal coefficients (DMDC) (see Figure 2). First, a model with a DMDC linear system is apt to produce relatively larger errors than that with a diagonally-dominant linear system, and is more sensitive to disturbances of the errors from the time-space discretization and the floating point operation. Second, if the simulated mixed flows are strongly dynamic (e.g., with frequent flow-regime conversions), the errors of a model with a DMDC linear system are possibly magnified, resulting in nonphysical oscillations of solutions at the transition points or even spurious flow-regime conversions.
The possible DMDC problem of the current linear solver is analyzed using three kinds of mixed flows which are differentiated with two characteristics. The first is “with flow-regime conversions (C1)” and the second is “in a closed pipe with wide-top cross-sections (C2)”.
TYPE I (without C1 or C2): Pressurized flows meet free-surface flows in nonclosed channels (manholes and open channels, e.g., cell i + 4 in Figure 2). In this case, although the DMDC may be significant, the regimes of the free-surface flow in nonclosed channels and the pressurized flow in closed pipes are both kept, without flow-regime conversions. TYPE II (with C1): Pressurized flows meet free-surface flows in a closed pipe with shrinking-top cross-sections (e.g., the transition point in Figure 1). Despite flow-regime conversions, the main-diagonal coefficients may be still assumed to be gradually varying around the transition points (the DMDC is not significant) if the unsteady property of the mixed flow is not strong. TYPE III (with C1 and C2): Pressurized flows meet free-surface flows in a closed pipe with wide-top cross-sections (e.g., the transition point in Figure 2), where flow-regime conversions and the DMDC are both significant.
In preliminary tests, the new model is revealed to be able to simulate the mixed flows of TYPE I and II stably. When the model is used to simulate TYPE-III mixed flows, it produces nonphysical oscillations or becomes unstable. The instability of the linear solver, arising from the simulation of the mixed flows with frequent flow-regime conversions in closed pipes with wide-top cross-sections, is called the DMDC problem. It is inferred that: in the nonlinear solver of [23], the nonlinear system has been continually updated using new flow variables in the iterative solutions and the potential instability (from similar discontinuous problems) is eliminated in a timely manner.
For simulating TYPE III mixed flows, an artificial slot (A-slot) is proposed to alleviate the possible DMDC problem of the linear solver and enhance its stability. Because the wet surface width of the cross-section is expressed explicitly in Equation (7), applying the A-slot to Equation (7) can be done straightforwardly:
C i a η i 1 n + 1 + C i b η i n + 1 + C i c η i + 1 n + 1 = M a x ( B i n , ε W i ) Δ x i η i n + r i n
where C i b = M a x ( B i n , ε W i ) Δ x i C i a C i c ; and ε is the percent that the width of the A-slot account for the representative width of the cross-section (Wi).
Equation (8) is coded instead of Equation (7) in the model, with the ε being an adjustable parameter according to the types of mixed flows. However, the DMDC is not obvious in simulations of many mixed flows (e.g., TYPE I and -II mixed flows), when the A-slot is not necessary and ε is set to 0. Hence, the A-slot is not an essential part of the current model. That is to say: the model without an A-slot can be applied to all type of mixed flows except the TYPE III mixed flows; the A-slot is only needed for the simulation of TYPE III mixed flows.

3. Results

The new model is demonstrated by tests of various mixed flows, including the aforementioned mixed flows of TYPE I, II and III. The simulation results using the new model are compared with analytical solutions, experiment data and results reported by other researchers.

3.1. Test Case 1

This test describes a pressurized flow in a horizontal pipe, connecting two reservoirs (see Figure 3). The pipe length (L) is 400 m and its square cross-section has an area of A = 1 m2. In the middle of the pipe, there is a valve that is closed. The water levels of the two reservoirs (η1 and η2) are kept constant, and their difference is η1η2 = 1 m. At time t = 0, the valve is opened suddenly and fully. Under the assumption that the pipe is frictionless, the water velocity and the pressure in the pipe can be described by analytical solutions as follows:
u ( x , t ) = u 0 tanh ( t / t 0 )
η ( x , t ) = η 2 + ( η 1 η 2 ) L x L cosh 2 ( t / t M )
where u 0 = 2 g ( η 1 η 2 ) , and t 0 = 2 L / u 0 . The analytical solutions are independent from shape, size, and depth of the pipe, as long as the pipe is fully filled. The u0 is calculated to be 4.43 m/s when η1η2 = 1 m.
The unsteady flow of this test falls into the category of TYPE-I mixed flows. The present model neglects water compressibility and pipe elasticity, which means the same discharge through the pipe from upstream to downstream. The model is tested on gradually reduced uniform grids (Meshes 1–5) whose grid scales are sequentially 40, 20, 16, 10, and 5 m, to establish the grid independence of solutions. In the simulations, Δt = 1 s and θ = 1, which are common in existing simulations. The model completes all the simulations stably. When the grid scale becomes smaller, the simulation result gradually converges at the analytical solution, as shown in Figure 4. The grid scale equal to or smaller than 16 m (Mesh 3) provides grid-independent results, when the errors in the simulated histories of the velocity and the pressure are all minor.

3.2. Test Case 2

The laboratory experiment in a closed pipe [32] is tested to demonstrate the model in simulating the free-surface and the pressurized flows connected by extreme forms (hydraulic jump). The pipes have a circular cross-section with an inner diameter of 0.22 m and consist of three reaches. The upstream and downstream reaches (L1 and L3) are both horizontal and 2 m long, which are connected by a downward inclined reach (L2) at an angle α = −10° having a length of 4 m. In an experiment, a constant water discharge of 0.03 m3/s and a constant pressure head of 0.554 m are imposed, respectively, at the upstream and downstream boundaries. The steady flow of the test falls into the category of TYPE-II mixed flows. A 1D grid of 0.01 m is used, and the implicit factor θ is set to 1.0. The wall friction in horizontal reaches of the pipe is neglected, which is the same as those in existing simulations [23,33]. Moreover, in the downward inclined reach, energy loss associated with air pockets [32,33] is assumed to act as an additional wall friction and be imposed on the free-surface parts, where the corresponding nm is determined by calibration.
In the first group of simulations (Δt = 0.003 s), to observe the evolution of mixed flows and the formation of the hydraulic jump, the discharge at the inlet of the pipe is gradually increased from 0 to 0.03 m3/s within a prescribed period (6 s), and is preserved at 0.03 m3/s after this period. The initial profile (at t = 0) and the simulated profiles (at t = 3, 6, 9, 12, 66 s) of the pressure are plotted in Figure 5. In the final steady state of the simulated mixed flow, the hydraulic jump occurs at about 0.3 m from the last measuring station. The results (at t ≥ 66 s) simulated by the current model agree well with the measured data and the simulation results reported by [23].
In the second group of simulations, the model is tested on gradually reduced time steps which are sequentially 0.01, 0.005, 0.002, 0.001, and 0.0005 s, to establish the time-step independence of the solutions. Correspondingly, the number of substeps for the ELM backtracking (Nbt) is, respectively, set to 20, 10, 4, 2, and 1. Using these time steps, a time-step sensitivity study of the current model is carried out, in terms of stability and accuracy.
On the one hand, the largest velocity of the simulations is found to be 3.28 m/s and occurs in the downward inclined reach. The CFL is calculated to be 3.38, 1.64, 0.66, 0.33, and 0.16, for the simulations using a Δt from 0.01 to 0.0005 s. The linear solver completes all the simulations stably. Moreover, it is found that the model begin to produce oscillated solutions when the Δt is equal to or larger than 0.1 s (CFL ≥ 3.38). In this test (a TYPE-II mixed flow), a gradually variation of the wet surface width at the circular cross-section can be observed in Figure 6. It means that the value of the main-diagonal coefficients of the linear system is gradually varying around the transition point (namely the DMDC problem is minor), so the A-slot is not needed. On the other hand, the time step, which is less than or equal to 0.005 s, provides time-step-independent results.

3.3. Test Case 3

The frictionless oscillating flows within a U-shaped tube are simulated. The tube has a square cross-section area of 1 m2. The initial conditions are mathematically defined by:
u ( x , 0 ) = 0
η ( x , 0 ) = z 0 + 0.01 cos ( π x / L )
where z0 is a reference level, and L is the horizontal length of the tube (L = 32 m), as shown in Figure 7. A 1D grid of Δx = 1 m is used and Δt is set to 0.01 s. Three scenarios of oscillating flows are considered below.
In Scenario 1, by setting z0 = 0.011, the oscillating flow in the horizontal section of the tube is fully pressurized and its oscillation frequency is given by 2 g / L . The flow of Scenario 1 falls into the category of TYPE-I mixed flows, and the A-slot is not needed. In Scenario 2, by setting z0 = –0.011, the oscillating flow is free-surface everywhere and its oscillation frequency is given by π g H / L , where H is the average water depth (0.989 m). The flow of Scenario 2 has only one flow regime and does not fall into any category of the TYPE-I, -II or -III mixed flows, and the A-slot is also not needed. In simulations of Scenarios 1 and 2, water-level histories (at x = 0) produced by the linear solver (LS) are plotted in Figure 8, and they are almost the same as the analytical solutions.
In Scenario 3, by setting z0 = 0, the horizontal section of the tube is filled with a partially pressurized and partially free-surface flow, which falls into the category of TYPE-III mixed flows. Moreover, the advection of the flow in this case is so weak that its solution provides little help to suppress the possible nonphysical oscillations in the solution of a linear system with the DMDC. As a result, the oscillating flow of Scenario 3 forms a challenging test case of TYPE-III mixed flows, to which the theoretical solutions are not available. The solutions of Scenario 3 have been only reported in studies by numerical methods, and the simulation results of Casulli and Stelling [23] is taken here as a reference. The new model is tested using six A-slots whose ε are sequentially 0.02, 0.03, 0.04, 0.05 0.1, and 0.2, to clarify its performance in simulating TYPE-III mixed flows.
For the six simulations of Scenario 3, the water-level histories (at x = 0) simulated by the linear solver (LS) using different A-slots are plotted in Figure 9 and compared with that in the reference [23]. The slots of ε < 0.05 are found to be unable to fully eliminate the nonphysical oscillations. When ε = 0.05–0.1, simulation results of the current model achieve the best agreement with those in [23]. When the ε is larger than 0.1, solutions of the linear solver begin to become diffusive. The suitable value of ε, which helps the current model achieve stable and accurate simulations of the TYPE-III mixed flows, is suggested to be 0.05–0.1.

4. Discussion

4.1. Differences between the New Model and Existing Models

In the introduction, two categories of models for the mixed flows have been reviewed, respectively, the TSE model (using two sets of equations) and the ASE model (using a single set of equations). The algorithms of the TSE models are often case-specific so that it is difficult to apply them to real applications. The ASE models can be further divided into two sub kinds, respectively, the ASE-1 (the pressure is often not explicitly included in the governing equations) and the ASE-2 (the pressure is often explicitly included). The new model solves two regimes of mixed flow using a unified formulation and the pressure term is explicitly included in the governing equations, so it is simpler than most of the existing TSE and ASE-1 models.
The new model can be classified to the ASE-2 models and may be considered as a variant of the ASE-2 model of Casulli and Stelling [23]. In the new model, the velocity–pressure coupling is solved by constructing and solving a linear system instead of the nonlinear system in [23]. The resulting linear system can be solved using a direct method and no iterations are needed. The new model is different with the model of [23], in terms of the continuity equation and its discretization (see Section 2.3.1), the construction of velocity-pressure coupling and the solution of algebraic equations (see Section 2.3.2). Generally, construction and solution of a nonlinear system is much more complex than those of a linear system, and the new model is obviously simpler than that of [23].
Moreover, the new model is revealed in tests to be applicable to various types of mixed flows, where stable and accurate simulations can be achieved at large time steps (CFL > 1). The ability of the new model is attractive for real applications. Moreover, the simplicity of the model allows that other researchers can develop a similar mixed-flow model easily and quickly.

4.2. The DMDC Problem and the A-Slot Approach

The main disadvantage of the new model is the possible DMDC problem, and this problem is clarified as follows: We analyze the characteristics of various typical mixed flows (see Section 2.3.3), and only the simulation of TYPE-III mixed flows (mixed flows with frequent conversions of flow regimes in a closed pipe with wide-top cross-sections) may trigger the DMDC problem of the linear solver. At the same time, the DMDC problem of the linear solver is not obvious in the simulations of other mixed flows (e.g., the TYPE I and -II mixed flows). The A-slot technique is proposed to alleviate the possible DMDC problem of the linear solver (eliminate the nonphysical oscillations of solutions), when the model is used to simulate the TYPE-III mixed flows.
On the one hand, the A-slot technique and the traditional Preissmann slot technique [20] can be differentiated in two aspects. First, these two slot techniques have quite different physical meanings. The A-slot technique is proposed to alleviate the possible instabilities of the linear solver related to the DMDC problem. The traditional Preissmann slot is used to help realizing “only one flow type (free-surface flow) throughout the whole pipe [9]” in the solution of mixed flows. Second, the model without an A-slot can be applied to all types of mixed flows except TYPE III mixed flows, and the A-slot is only needed for the simulation of TYPE III mixed flows. As a result, the A-slot is not an essential part of the new model except that a special case (the TYPE-III mixed flows) is simulated. To be different, the Preissmann slot technique is an essential part of the model no matter which kind of mixed flows is simulated.
On the other hand, an interesting finding is that the A-slot and the Preissmann slot favor quite similar slot widths. Although the slot width (ε) of a Preissmann slot should be sized to match the assumed celerity in the pressurized part of the flow, in most cases the ε is arbitrarily chosen based on numerical considerations [17]. In existing studies, the ε of a Preissmann slot is often set to 0.02 [13,18], 0.05 [2,34], 0.1 [14,16], etc. In tests of this study, the suitable value of ε, which enables the linear solver to achieve stable and accurate simulations of TYPE-III mixed flows, is suggested to be 0.05–0.1, and they are quite similar to the widths of the Preissmann slots used in existing models.

5. Conclusions

A semi–implicit numerical model with a linear solver is proposed for the free-surface and pressurized mixed flows in hydraulic systems. It solves the two flow regimes within a unified formulation and is much simpler than existing similar models for mixed flows. Using a local linearization and an ELM, the new model only needs to solve a tridiagonal linear system arising from the velocity-pressure coupling. The linear system has a symmetric and positive definite coefficient matrix, and can be directly solved without iterations.
Characteristics of three representative types (TYPE I, II, and III) of mixed flows are analyzed, and the applicability of the new model to each type of them is discussed. TYPE-III mixed flow, characterized by both flow-regime conversions and a closed pipe with wide-top cross-sections, is the most challenging case. In simulations of TYPE-III mixed flows using the linear solver, an artificial slot (A-slot) technique is proposed to cope with the possible instabilities related to the discontinuous main-diagonal coefficients of the linear system.
The new model is tested using various mixed flows of TYPE I, II, and III, where the simulation results agree with the analytical solutions, experiment data, and the simulation results reported by former researchers. Sensitivity studies of grid scales and time steps are both performed, where a common grid scale provides grid-independent results and a common time step provides time-step-independent results. Moreover, the model is revealed to achieve stable and accurate simulations at large time steps for which the CFL can be greater than 1. In simulations of a TYPE-III mixed flow using an A-slot, a slot-width sensitivity study is also performed. To ensure stable and accurate simulations of TYPE-III mixed flows, the suitable slot-width ratio (ε) for the linear solver is suggested to be 0.05–0.1.

Author Contributions

Conceptualization: D.H. and Z.J.; methodology: D.H. and Z.J.; software: D.H. and S.L.; validation: D.H., S.L. and Z.J.; formal analysis: D.H., S.L. and Z.J.; investigation: D.H. and Z.J.; resources: S.Y., and Z.J.; data curation: S.Y. and D.H.; writing—original draft preparation: D.H. and Z.J.; writing—review and editing: S.L. and D.H.; visualization: S.L.; supervision: S.Y. and Z.J.; project administration: S.Y. and Z.J.; funding acquisition: D.H., S.Y. and Z.J.

Funding

Financial support from Science and Technology Program of Guizhou Province, China (grant no. 2017-3005-4), the Fundamental Research Funds for the Central Universities (2017KFYXJJ197), China’s National Natural Science Foundation (51339001, 51779015, and 51479009) are acknowledged.

Conflicts of Interest

The authors declare no conflicts of interest.

Notation

The following symbols and abbreviations are used in this paper:
AWetted cross-sectional area (conveyance area);
ApThe area of the cross-section;
A-slotArtificial slot for coping with the discontinuous main-diagonal coefficients problem
B= A / η , the wet surface width of a cross-section for the free-surface parts and is equal to 0 for the pressurized parts of a closed pipe;
gGravity acceleration;
neThe quantities of cells;
nsThe quantities of sides;
nmManning’s roughness coefficient;
NbtThe number of substeps of the ELM
RHydraulic radius;
tTime;
uAveraged velocity of cross-section;
ubtSolution of the advection term uisng the ELM;
WiRepresentative width of the cross-section;
xLongitudinal distance along the channel;
zcThe celling of the cross-section;
tTime;
ΔtTime step
ΔxiThe length of cell i;
Δxi+1/2The distance between the centers of cells i and i + 1;
θImplicit factor
ηWater level measured from an undisturbed reference water surface (for free-surface flows); piezometric head measured from an undisturbed reference water surface (for pressurized flows);
εThe percent that the width of the A-slot account for the representative width
ASEA single set of equations;
CLOA mark which is used to distinguish closed and nonclosed reaches;
DMDCDiscontinuous main-diagonal coefficients
ELMEulerian–Lagrangian method;
PREA mark which is used to distinguish free-surface and pressurized reaches;
SVESaint–Venant equations
TSETwo sets of equations;

References

  1. Yen, B.C. Hydraulics of sewer systems. In Stormwater Collection Systems Design Handbook; McGraw-Hill: New York, NY, USA, 2001. [Google Scholar]
  2. Kerger, F.; Archambeau, P.; Erpicum, S.; Dewals, B.J.; Pirotton, M. A fast universal solver for 1D continuous and discontinuous steady flows in rivers and pipes. Int. J. Numer. Meth. Fluids 2011, 66, 38–48. [Google Scholar] [CrossRef]
  3. Hamam, M.A.; McCorquodale, J.A. Transient conditions in the transition from gravity to surcharged sewer flow. Can. J. Civ. Eng. 1982, 9, 189–196. [Google Scholar] [CrossRef]
  4. Li, J.; McCorquodale, A. Modeling mixed flow in storm sewers. J. Hydraul. Eng. 1999, 125, 1170–1180. [Google Scholar] [CrossRef]
  5. Zhou, F.; Hicks, F.E.; Steffler, P.M. Transient flow in a rapidly filling horizontal pipe containing trapped air. J. Hydraul. Eng. 2002, 28, 25–634. [Google Scholar] [CrossRef]
  6. Song, C.S.S.; Cardle, J.A.; Leung, K.S. Transient mixed flow models for storm sewers. J. Hydraul. Eng. 1983, 109, 1487–1504. [Google Scholar] [CrossRef]
  7. Cardle, J.A.; Song, C.S.S. Mathematical modeling of unsteady flow in storm sewers. Int. J. Eng. Fluid Mech. 1988, 1, 495–518. [Google Scholar]
  8. Guo, Q.; Song, C.S.S. Surging in urban storm drainage systems. J. Hydraul. Eng. 1990, 116, 1523–1537. [Google Scholar] [CrossRef]
  9. Fuamba, M. Contribution on transient flow modelling in storm sewers. J. Hydraul. Res. 2002, 40, 685–693. [Google Scholar] [CrossRef]
  10. Wang, K.H.; Shen, Q.; Zhang, B. Modeling propagation of pressure surges with the formation of an air pocket in pipelines. Comput. Fluids 2003, 32, 1179–1194. [Google Scholar] [CrossRef]
  11. Politano, M.; Odgaard, A.J.; Klecan, W. Case study: Numerical evaluation of hydraulic transients in a combined sewer overflow tunnel system. J. Hydraul. Eng. 2007, 133, 1103–1110. [Google Scholar] [CrossRef]
  12. Bousso, S.; Daynou, M.; Fuamba, M. Numerical Modeling of Mixed Flows in Storm Water Systems: Critical Review of Literature. J. Hydraul. Eng. 2013, 139, 385–396. [Google Scholar] [CrossRef]
  13. Garcia-Navarro, P.; Priestley, A.; Alcrudo, S. Implicit method for water flow modeling in channels and pipes. J. Hydraul. Res. 1994, 32, 721–742. [Google Scholar] [CrossRef]
  14. Capart, H.; Sillen, X.; Zech, Y. Numerical and experimental water transients in sewer pipes. J. Hydraul. Res. 1997, 35, 659–672. [Google Scholar] [CrossRef]
  15. Ji, Z. General hydrodynamic model for sewer/channel network systems. J. Hydraul. Eng. 1998, 124, 307–315. [Google Scholar] [CrossRef]
  16. Trajkovic, B.; Ivetic, M.; Calomino, F.; D’Ippolito, A. Investigation of transition from free surface to pressurized flow in a circular pipe. Water Sci. Technol. 1999, 39, 105–112. [Google Scholar] [CrossRef]
  17. Vasconcelos, J.G.; Wright, S.J.; Roe, P.L. Improved Simulation of Flow Regime Transition in Sewers: Two-Component Pressure Approach. J. Hydraul. Eng. 2006, 132, 553–562. [Google Scholar] [CrossRef]
  18. León, A.S.; Ghidaoui, M.S.; Schmidt, A.R.; García, M.H. Application of Godunov-type schemes to transient mixed flows. J. Hydraul. Res. 2009, 47, 147–156. [Google Scholar] [CrossRef]
  19. Bourdarias, C.; Gerbi, S. A conservative model for unsteady flows in deformable closed pipes and its implicit second-order finite volume discretisation. Comput. Fluids 2008, 37, 1225–1237. [Google Scholar] [CrossRef]
  20. Cunge, J.A.; Wegner, M. Numerical integration of Barre de Saint-Venant’s flow equations by means of implicit scheme of finite differences. Houille Blanche 1964, 19, 33–39. [Google Scholar] [CrossRef]
  21. Sanders, B.F.; Bradford, S.F. Network implementation of the two-component pressure approach for transient flow in storm sewers. J. Hydraul. Eng. 2011, 137, 158–172. [Google Scholar] [CrossRef]
  22. Aldrighetti, E. Computational Hydraulic Techniques for the Saint Venant Equations in Arbitrarily Shaped Geometry. Ph.D. Thesis, University of Trento, Trento, Italy, 2007. [Google Scholar]
  23. Casulli, V.; Stelling, G.S. A semi-implicit numerical model for urban drainage systems. Int. J. Numer. Meth. Fluids 2013, 73, 600–614. [Google Scholar] [CrossRef]
  24. Dumbser, M.; Iben, U.; Ioriattia, M. An efficient semi-implicit finite volume method for axially symmetric compressible flows in compliant tubes. Appl. Numer. Math. 2015, 89, 24–44. [Google Scholar] [CrossRef]
  25. Leskens, J.G.; Brugnach, M.; Hoekstra, A.Y.; Schuurmans, W. Why are decisions in flood disaster management so poorly supported by information from flood models? Environ. Model. Softw. 2014, 53, 53–61. [Google Scholar] [CrossRef]
  26. Storm Water Management Model (SWMM). Storm Water Management Model Reference Manual Volume II-Hydraulics; EPA/600/R-17/111; Rossman: Los Angeles, CA, USA; National Risk Management Laboratory: Cincinnati, OH, USA, 2017.
  27. Dimou, K. 3-D Hybrid Eulerian–Lagrangian/Particle Tracking Model for Simulating Mass Transport in Coastal Water Bodies. Ph.D. Thesis, Department of Civil Engineering, Massachusetts Institute of Technology, Cambridge, MA, USA, 1992. [Google Scholar]
  28. Casulli, V.; Zanolli, P. Semi-implicit numerical modeling of nonhydrostatic free-surface flows for environmental problems. Math. Comput. Model. 2002, 36, 1131–1149. [Google Scholar] [CrossRef]
  29. Zhang, Y.L.; Baptista, A.M.; Myers, E.P. A cross-scale model for 3D baroclinic circulation in estuary-plume-shelf systems, I. Formulation and skill assessment. Cont. Shelf Res. 2004, 24, 2187–2214. [Google Scholar] [CrossRef]
  30. Hu, D.C.; Zhong, D.Y.; Zhang, H.W.; Wang, G.Q. Prediction–Correction Method for Parallelizing Implicit 2D Hydrodynamic Models I Scheme. J. Hydraul. Eng. 2015, 141, 04015014. [Google Scholar] [CrossRef]
  31. Hu, D.C.; Zhong, D.Y.; Zhu, Y.H.; Wang, G.Q. Prediction–Correction Method for Parallelizing Implicit 2D Hydrodynamic Models II Application. J. Hydraul. Eng. 2015, 141, 06015008. [Google Scholar] [CrossRef]
  32. Pothof, I. Co–Current Air–Water Flow in Downward Sloping Pipes; Gildeprint Drukkerijen BV: Enschede, The Netherland, 2011; ISBN 978-90-8957-018-5. [Google Scholar]
  33. Lubbers, C.L. On Gas pockets in Wastewater Pressure Mains and Their Effect on Hydraulic Performance. Ph.D. Thesis, Department of Civil Engineering and Geosciences, University of Technology, Delft, The Netherlands, 2007. [Google Scholar]
  34. Ferreri, G.B.; Freni, G.; Tomaselli, P. Ability of Preissmann slot scheme to simulate smooth pressurisation transient in sewers. Water Sci. Technol. 2010, 62, 1848–1858. [Google Scholar] [CrossRef]
Figure 1. Sketch of mixed flow in a closed sewer pipe with manholes and computational grids.
Figure 1. Sketch of mixed flow in a closed sewer pipe with manholes and computational grids.
Water 11 01979 g001
Figure 2. Description of discontinuous distribution of the main-diagonal coefficients of linear system (using the case of a closed pipe with rectangular cross-sections).
Figure 2. Description of discontinuous distribution of the main-diagonal coefficients of linear system (using the case of a closed pipe with rectangular cross-sections).
Water 11 01979 g002
Figure 3. Two reservoirs connected by a horizontal pressurized pipe.
Figure 3. Two reservoirs connected by a horizontal pressurized pipe.
Water 11 01979 g003
Figure 4. Simulated histories at the first grid point inside the pipe (using different grid scales): (a) the velocity; (b) the pressure head.
Figure 4. Simulated histories at the first grid point inside the pipe (using different grid scales): (a) the velocity; (b) the pressure head.
Water 11 01979 g004
Figure 5. Evolution of mixed flows in closed circular pipe and formation of hydraulic jump.
Figure 5. Evolution of mixed flows in closed circular pipe and formation of hydraulic jump.
Water 11 01979 g005
Figure 6. Simulated profiles of pressure head and hydraulic jump using different Δt.
Figure 6. Simulated profiles of pressure head and hydraulic jump using different Δt.
Water 11 01979 g006
Figure 7. U-shaped tube with a square cross-section.
Figure 7. U-shaped tube with a square cross-section.
Water 11 01979 g007
Figure 8. Simulated histories of water levels in scenarios of pure pressurized/free-surface flows.
Figure 8. Simulated histories of water levels in scenarios of pure pressurized/free-surface flows.
Water 11 01979 g008
Figure 9. Simulated histories of water levels in scenarios of mixed flows with an A-slot.
Figure 9. Simulated histories of water levels in scenarios of mixed flows with an A-slot.
Water 11 01979 g009

Share and Cite

MDPI and ACS Style

Hu, D.; Li, S.; Yao, S.; Jin, Z. A Simple and Unified Linear Solver for Free-Surface and Pressurized Mixed Flows in Hydraulic Systems. Water 2019, 11, 1979. https://doi.org/10.3390/w11101979

AMA Style

Hu D, Li S, Yao S, Jin Z. A Simple and Unified Linear Solver for Free-Surface and Pressurized Mixed Flows in Hydraulic Systems. Water. 2019; 11(10):1979. https://doi.org/10.3390/w11101979

Chicago/Turabian Style

Hu, Dechao, Songping Li, Shiming Yao, and Zhongwu Jin. 2019. "A Simple and Unified Linear Solver for Free-Surface and Pressurized Mixed Flows in Hydraulic Systems" Water 11, no. 10: 1979. https://doi.org/10.3390/w11101979

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop