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

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.


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 Water 2019, 11,1979; doi:10.3390/w11101979 www.mdpi.com/journal/water Water 2019, 11,1979 2 of 15 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.

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

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: 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; n m = 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: 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 (A p ) when the water level is equal to or larger than the celling of the cross-section (z c ). 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.

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 Water 2019, 11,1979 4 of 15 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, ∆x i represents the length of the cell; B i is the wet surface width and is related to the cell water level that varies with respect to time; A i 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, ∆x i+1/2 is the distance between the centers of cells i and i + 1; A i+1/2 and R i+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. 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.

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 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

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 bt , where "N bt " 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 u bt . The cell-face velocity is then directly updated with u bt , 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): where θ is the implicit factor; ∆t is the time step; superscript n indicates time level n; and the riverbed friction term is discretized using u bt and u n+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: ; G n i+1/2 is the incorporated explicit terms and G n i+1/2 /F n i+1/2 can be regarded as the provisional solution by solving all the explicit terms, 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 n+1 i and the A n i 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 η n+1 i and η n i 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 ∂(Au)/∂x accounts for the calculation of the water fluxes through cell faces. For the discretization of the ∂(Au)/∂x, we can use either the A n i±1/2 or the A n+1 i±1/2 (at cell faces), which may only influence the order of the accuracy of solutions. In practice, the A n i±1/2 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 (B i ) is explicitly represented by the value of time level n (B n i ) 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: The nonlinear terms (B i and A i±1/2 ) are present, respectively, on the left hand side (LHS) and the right hand side (RHS) of Equation (6). However, when the B i on the LHS and the A i±1/2 on the RHS have been explicitly represented by B n i and A n i±1/2 (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 η n+1 i − A i η n i , 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.

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): (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 a i < 0, C c i < 0 and C b i > 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 u n+1 defined at a side center is calculated by Equation (5), and then the A i and the R i 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.

Considerations of Abrupt Change of Wet Surface Width in Mixed Flows
The wet surface width of a cross-section (B i ) has a positive finite value for free-surface flows, and is equal to 0 for pressurized flows in pipes. Hence, abrupt changes of the B i 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 B i changes abruptly from the width of W i (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 B i 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 B i 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 B i has a positive value and the coefficient matrix of the linear system arising from the velocity-pressure coupling is diagonally dominant. The larger the B i is, the larger the main-diagonal coefficient is. When the model is used to simulate pure pressurized flows, with B i = 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 B i 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.
Water 2019, 11, x FOR PEER REVIEW 7 of 16 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 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: where C b i = Max(B n i , εW i )∆x i − C a i − C c i ; and ε is the percent that the width of the A-slot account for the representative width of the cross-section (W i ).
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.

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.

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 m 2 . 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: where u 0 = 2g(η 1 − η 2 ), and t 0 = 2L/u 0 . The analytical solutions are independent from shape, size, and depth of the pipe, as long as the pipe is fully filled. The u 0 is calculated to be 4.43 m/s when η 1 − η 2 = 1 m. width of the cross-section is expressed explicitly in Equation (7), applying the A-slot to Equation (7) can be done straightforwardly: 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.

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.

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 m 2 . 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:  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 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.

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 m 3 /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 m 3 /s within a prescribed period (6 s), and is preserved at 0.03 m 3 /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].

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 m 3 /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 n m 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 m 3 /s within a prescribed period (6 s), and is preserved at 0.03 m 3 /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 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 (N bt ) 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.

Test Case 3
The frictionless oscillating flows within a U-shaped tube are simulated. The tube has a square cross-section area of 1 m 2 . The initial conditions are mathematically defined by: 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 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.

Test Case 3
The frictionless oscillating flows within a U-shaped tube are simulated. The tube has a square cross-section area of 1 m 2 . The initial conditions are mathematically defined by: η(x, 0) = z 0 + 0.01 cos(πx/L) (12) where z 0 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 z 0 = 0.011, the oscillating flow in the horizontal section of the tube is fully pressurized and its oscillation frequency is given by 2g/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 z 0 = -0.011, the oscillating flow is free-surface everywhere and its oscillation frequency is given by π gH/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 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 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. In Scenario 3, by setting z 0 = 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. 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. x(m) 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.

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.

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.

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.

Conflicts of Interest:
The authors declare no conflicts of interest.

Notation
The following symbols and abbreviations are used in this paper:

A
Wetted cross-sectional area (conveyance area); A p The area of the cross-section; A-slot Artificial 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; g Gravity acceleration; ne The quantities of cells; ns The quantities of sides; n m Manning's roughness coefficient; N bt The number of substeps of the ELM R Hydraulic radius; t Time; u Averaged velocity of cross-section; u bt Solution of the advection term uisng the ELM; W i Representative width of the cross-section; x Longitudinal distance along the channel; z c The celling of the cross-section; t Time; ∆t Time step ∆x i The length of cell i; ∆x i+1/2 The 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 ASE A single set of equations; CLO A mark which is used to distinguish closed and nonclosed reaches; DMDC Discontinuous main-diagonal coefficients ELM Eulerian-Lagrangian method; PRE A mark which is used to distinguish free-surface and pressurized reaches; SVE Saint-Venant equations TSE Two sets of equations;