1D-3D coupled simulation method of hydraulic transients in ultra-long hydraulic systems based on OpenFOAM

Abstract

Transients in hydraulic systems are normally simulated by one-dimensional (1D) method of characteristics (MOC). When local three-dimensional (3D) flow patterns have crucial impacts on the hydraulic system characteristics, 3D computational fluid dynamics (CFD) simulation should be incorporated. Since simulating an entire long system by the 3D CFD method is unrealistic, we have proposed several 1D-3D coupled simulation methods. However, the coupled methods are difficult to apply and popularise because most CFD codes are inaccessible and implementation details are unclear. Consequently, we utilised the open source toolbox OpenFOAM to develop a 1D-3D coupled simulation method and proved its feasibility by applying it to analysing practical problems. This paper describes the implementations of water compressibility, 1D-3D coupling, gas-liquid flow simulation, and dynamic mesh; shows their verifications of accuracy by benchmark cases; demonstrates the simulations of the closing and opening processes of a 194 km long water conveyance system. The gas-liquid flow patterns in the gate chamber and the maximum/minimum water head envelopes along the tunnels, which support the design optimisation of the project, are presented. It is shown that the method is computationally efficient and accurate, and can be applied to hydraulic transient analyses for water conveyance systems, pumped-storage power stations, etc.

KEYWORDS:

• 1D-3D coupling
• MOC
• CFD
• hydraulic transients
• gas-liquid flow
• long hydraulic system

1. Introduction

Nowadays, many water conveyance systems, hydropower stations, pumping stations, pumped-storage power stations, and other hydraulic projects are under construction in the world. Ensuring safety and minimising investment are the main concerns of design and optimisation. Because hydraulic transients are related to safety, the study of transient processes is indispensable for most hydraulic systems. Some local flow fields in hydraulic systems, such as those in pump-turbines of pumped-storage systems and control gates of water conveyance systems, have important impacts on the operating characteristics of the whole systems. The histories over time and distributions along pipelines of discharge and head are closely related to the local flow patterns. Therefore, the coupled simulations of the local flow fields and water hammer in pipelines should be performed in the hydraulic transient study of such systems.

For decades, the main approach to studying the hydraulic transient processes is one-dimensional (1D) numerical simulation, where using the method of characteristics (MOC) is dominant. However, the 1D methods cannot consider the three-dimensional (3D) flow patterns and their impact on the hydraulic system characteristics. When the interaction between local flow and system is significant, multi-dimensional computational fluid dynamics (CFD) should be used. 3D CFD method has been used for simulating water hammer problems with turbulence flow, and its accuracy has been validated by comparing with experiments (Warda et al., 2020; Yang et al., 2017). Muti-dimensional water hammer phenomena can also be observed (Zhang et al., 2022). The local flow characteristics of multiphase flow (He et al., 2022; Ma et al., 2020; Qi et al., 2022; Xia et al., 2013), cavitation (Tang et al., 2020; Warda et al., 2020), moving boundary (Li et al. 2019; Ye et al., 2022), etc., can be well provided. However, in practical engineering applications, the scale of the simulation object is often very large, for example, an inter-basin water conveyance system is generally tens or even hundreds of kilometres long (Nan et al., 2018; Zhu et al., 2012). In such cases, simulating the whole system by the 3D CFD method is unrealistic because of the unacceptable computational resources and time. Also, we found that using 3D CFD to evaluate the frictional head loss of long pipes is not accurate (Cai et al., 2016), but the head loss for long hydraulic systems is normally an important influence factor. The conclusion is that for systems with obvious 1D characteristics, 1D MOC simulation is competent (Liu et al., 2021; Urbanowicz et al., 2021), and full 3D CFD simulation is neither realistic nor necessary.

Therefore, to consider the interaction of local flow patterns with long pipeline systems, 1D-3D coupled simulation methods have been proposed, in which the long pressurised pipeline system is simulated by the 1D MOC and the local flow field is simulated by the 3D CFD method. The core technique is solving the boundary conditions on the coupling interface with the flow field information of both 1D and 3D sides simultaneously. There exist different kinds of coupling algorithms. Some used the flow information of the previous time instead of the current time (Ruprecht & Helmrich, 2003) to specify the boundaries of both sides. Because of its approximation, this method requires a very small time step. Some adopted iterative calculation (Huang et al., 2012; X. Li, Tang, et al., 2019; Wu et al., 2015), which leads to an increase in computation and difficulty in implementation. There are also explicit calculation methods, which are popular in recent years because of their high accuracy and efficiency, and easy implementation. For example, Zhang and Cheng (2012) proposed the partly overlapped coupling method, in which the 1D and 3D meshes are partly overlapped at the junction and a small segment of the 3D mesh is used instead of the 1D mesh, thus realising explicit calculation. Galindo et al. (2011) and Wang et al. (2017) proposed coupling methods in explicit format based on the Riemann invariant equation, in which the characteristic equation on the 1D side and the Riemannian invariant equation on the 3D side are used in the interface calculation. These 1D-3D coupled algorithms were used in the simulations of transient processes in practical engineering projects, such as hydropower stations (Mandair et al., 2021), pumping stations (Li et al. 2019) and pumped-storage power stations (Fu et al., 2023; Yang et al., 2020).

Most of the existing implementations of 1D-3D coupling are based on user-defined function (UDF) in commercial CFD software. The core codes of commercial software are inaccessible, leading to complicated secondary development and inconvenient popularisation. Moreover, the numerical models required for the coupled method and their implementation details are still not clearly described in existing literature, making it difficult for new researchers in this field to use the method.

To provide an open-source and easy-to-implement 1D-3D coupled method for efficient and accurate simulations of long hydraulic systems, we chose OpenFOAM (v2112) to develop the computing programme and showed the details of implementation. Then the method was applied to analyse the opening and closing processes of a practical long-distance water conveyance system. This paper is structured as follows. Section 2 describes the numerical models related to the 1D-3D coupled simulation and their implementation methods. Section 3 verifies the accuracy of the related algorithms by benchmark cases. In Section 4, the 1D-3D coupled algorithm is used to simulate the transient processes of a 194 km long water conveyance system, whose control gates are in the middle-lower reach.

2. Numerical models and implementation

This section describes the numerical models of the proposed 1D-3D coupled simulation method with their implementation for a practical long-distance water conveyance system. The basic numerical models include 1D MOC for pipelines, water compressibility for 3D CFD, and 1D-3D coupling for interface; the additional models include gas-liquid flow and dynamic mesh, which are specific for the practical example.

2.1 Basic numerical models for 1D-3D coupled simulation

2.1.1 1D method of characteristics (MOC)

The transient flow in a pressurised pipe can be simply described by the following continuity equation and momentum equation (Chaudhry, 2014). (1) $\begin{array}{rl}\frac{\mathrm{\partial }H}{\mathrm{\partial }t}+\frac{a^2}{gA}\frac{\mathrm{\partial }Q}{\mathrm{\partial }x}& =0\end{array}$(1) (2) $\begin{array}{rl}\frac{\mathrm{\partial }Q}{\mathrm{\partial }\text{t}}+gA\frac{\mathrm{\partial }H}{\mathrm{\partial }x}+\frac{fQ|Q|}{2DA}& =0\end{array}$(2) where x and t are the space and time coordinates, respectively; H is the piezometric head, and Q is the volumetric flow rate; a is the wave speed in water, g is the gravitational acceleration, and f is the Darcy friction loss coefficient; A and D are the cross-sectional area and diameter of the pipe, respectively.

Equations (1) and (2) form a quasilinear hyperbolic system of partial differential equations. To solve it, the MOC is used to convert it into a system of ordinary differential equations, which are then integrated along the positive and negative characteristic lines (Figure 1) to obtain the positive (Equation 3) and negative (Equation 4) characteristic equations (Chaudhry, 2014). (3) $\begin{array}{rl}{\mathrm{C}}^{+}:Q_{\mathrm{P}}& =C_p-C_a{H}_{\mathrm{P}}\end{array}$(3) (4) $\begin{array}{rl}{\mathrm{C}}^{-}:Q_{\mathrm{P}}& =C_n+C_a{H}_{\mathrm{P}}\end{array}$(4) where C_a is the pipe characteristic constant, expressed as C_a = gA / a; C_p and C_n are the positive and negative characteristic coefficients, respectively, with expressions C_p = Q_A + C_aH_A − ΔtQ_A |Q_A| f / 2DA and C_n = Q_B − C_aH_B − ΔtQ_B |Q_B| f / 2DA, which can be obtained from the information in the previous time step; the subscripts P, A and B denote the current computing node, the upstream node of the previous time step and the downstream node of the previous time step, respectively.

Figure 1. Schematic of 1D MOC.

So far, we can use Equations (3) and (4), combined with the pipe’s initial conditions, to calculate the discharge and head at each moment of the internal nodes of the pipe. Then, with the boundary conditions at the inlet and outlet of the pipe, the discharge and head at the boundary nodes of the pipe can be calculated. The boundary conditions depend on the type of pipe junctions, such as reservoirs, surge shafts and pipe connectors. In the case of a reservoir, the head is a fixed value; in the case of a surge shaft or pipe connector, it is required to solve the mass and energy conservation equations, and the Newton-Raphson iterative method is usually recommended.

Since 1D and 3D computations need to be performed in the same programme, and the 3D computing software OpenFOAM used in this study is programmed in C++, the 1D computing code should also be based on C++. In addition, for the sake of the universality and extensibility of the programme, it is recommended to adopt the object-oriented programming idea and define pipes and junctions as classes.

2.1.2 Water compressibility for 3D CFD

In the 1D-3D coupled simulation, the water compressibility that is related to wave speed in the 3D CFD domain must be considered, namely, compressible flow models should be used instead of incompressible ones where the wave speed is infinite.

The water compressibility should be set according to the wave speed in the 3D domain, which is consistent with that in the adjacent 1D domain. For water, the density and pressure are approximately linearly related; thus the linear equation of state can be used (Jablonská, 2014; Wang et al., 2017). (5) $\rho ={\rho }_0+\psi (p-p_0)$(5) where ρ is the current density; ρ₀ is the reference density, taken as 1000 kg m⁻³ in this study; p is the current pressure (relative pressure); p₀ is the reference pressure, taken as 0 Pa in this study; ψ is the compressibility, defined as ψ = dρ / dp. According to the bulk modulus of elasticity of water (K), we can obtain (6) $K=-\frac{\mathrm{d}p}{\mathrm{d}V/V}=\frac{\mathrm{d}p}{\mathrm{d}\rho /\rho }=\frac{\rho }{\psi }$(6) where V is the volume of water. The wave speed (a) satisfies the following equation (Ghidaoui et al., 2005) (7) $\frac{1}{a^2}=\frac{\rho }{K}+\frac{\rho }{Ee/D}\text{or}{a}^2=\frac{K/\rho }{1+KD/Ee}$(7) where E is the Young’s modulus of elasticity of the pipe material, e is the thickness of the pipe wall, and D is the Diameter of the pipe. Equation (7) indicates that the wave speed is affected by the water compressibility and pipe flexibility. However, the pipe wall mesh in 3D CFD simulations is fixed, i.e. rigid (E = ∞). Therefore, the wave speed in the 3D domain can only be reflected by the water compressibility, and its expression becomes (8) $\frac{1}{a^2}=\frac{\rho }{K}=\psi$(8) By substituting the value of a in the adjacent 1D domain into Equation (8), we can calculate ψ and define the equation of state (Equation (5)). In this study, given the pipe flexibility, the wave speed (a) is taken as 1200 m s⁻¹ in Section 4, which is slightly less than the speed of sound in water (about 1500 m s⁻¹).

In OpenFOAM, the definition of the equation of state is located in the ‘specie’ library of the ‘thermophsicalModels’ module. Equation (5) can be implemented with simple modifications based on the ‘perfectGas’ class template, focusing on modifying the ‘rho’ function. Note that the compressible flow solvers in OpenFOAM solve the energy equation by default, however, for the water hammer phenomenon in pipelines, which is the focus of this study, the effect of temperature is insignificant, and therefore, the solution of the energy equation can be blocked to speed up the computation. The specific approach to block the solution of the energy equation, taking the ‘rhoPimpleFoam’ solver as an example, is as follows: comment out most of the code in the file ‘EEqn.H’, and keep only the code ‘thermo.correct()’, which is used to update parameters such as density, making the equation of state work.

2.1.3. 1D-3D coupling

The 1D-3D coupling is essentially exchanges of data between the 1D and 3D computation. The calculation of the boundary conditions at the coupling interface requires contributions from both sides, i.e. the characteristic equation from the 1D side and the Riemann invariant equation from the 3D side (Galindo et al., 2011; Wang et al., 2017). The brief principles are as follows.

Neglecting the minor friction term fQ |Q| / 2DA in Equation (2), Equations (1) and (2) can be written in the following matrix form (Zhao & Ghidaoui, 2004), which can be solved as the Riemann problem. (9) $\frac{\mathrm{\partial }\mathit{u}}{\mathrm{\partial }t}+\frac{\mathrm{\partial }\mathit{f}(\mathit{u})}{\mathrm{\partial }x}=0,$(9) where $\mathit{f}(\mathit{u})=\mathrm{A}\mathit{u}$, with $\mathrm{A}=\left(\begin{array}{cc}0& a^2/gA\\ gA& 0\end{array}\right)$ and $\mathit{u}=\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}H\\ Q\end{array}\right)$. The eigenvalues and the corresponding right eigenvectors of matrix A are (10) ${\lambda }_1=a,{\mathit{p}}_1=({p_1}^{(1)},p_1^{(2)}{)}^{\mathrm{T}}=(a,gA{)}^{\mathrm{T}}$(10) (11) ${\lambda }_2=-a,{\mathit{p}}_2=(p_2^{(1)},p_2^{(2)}{)}^{\mathrm{T}}=(a,-gA{)}^{\mathrm{T}}$(11) According to the generalised Riemann invariants (Toro, 2013), on each of the two characteristic lines, we can obtain (12) $\frac{\text{d}{u}_1}{p_1^{(1)}}=\frac{\mathrm{d}{u}_2}{p_1^{(2)}}\Rightarrow \frac{\text{d}H}{a}=\frac{\mathrm{d}Q}{gA}$(12) (13) $\frac{\text{d}{u}_1}{p_2^{(1)}}=\frac{\mathrm{d}{u}_2}{p_2^{(2)}}\Rightarrow \frac{\mathrm{d}H}{a}=-\frac{\mathrm{d}Q}{gA}$(13) Integrating both sides of Equations (12) and (13), the Riemann invariant equations associated with the C⁺ and C⁻ characteristic lines can be obtained, respectively. (14) $\mathrm{R}{\mathrm{I}}^{-}:H-\frac{a}{gA}Q=\text{const}$(14) (15) $\mathrm{R}{\mathrm{I}}^{+}:H+\frac{a}{gA}Q=\mathrm{const}$(15) The contribution of the 3D side can be calculated with Equations (14) and (15), and then combined with the 1D side contribution obtained by the characteristic equation, the flow field information at the 1D-3D coupling interface can be obtained.

For a 1D-3D coupled example as shown in Figure 2, the explicit calculation procedures of the flow field information at the interface are as follows.

1. Specify the initial boundary conditions on the 3D side, and compute the 3D flow field at time t₀.

2. Extract the pressure and velocity at the volume centre of the 3D boundary layer mesh and take the average value. Then, convert them into the piezometric head and volumetric flow rate respectively, i.e. $H_{0.5}^{t_0}$ and $Q_{0.5}^{t_0}$.

3. According to the Riemann invariant equation RI⁻ (Equation (14)) associated with the negative characteristic line, obtain (16) $H_0^{t_1}-\frac{a}{gA}{Q}_0^{t_1}=H_{0.5}^{t_0}-\frac{a}{gA}{Q}_{0.5}^{t_0}$(16)

4. Specify the initial boundary conditions on the 1D side, and compute the 1D flow field at time t₀.

5. Extract the piezometric head and volumetric flow rate at the computing node adjacent to the 1D boundary, i.e. $H_{n-1}^{t_0}$ and $Q_{n-1}^{t_0}$.

6. According to the positive characteristic equation C⁺ (Equation (3)), obtain (17) $Q_n^{t_1}=C_p^{t_0}-C_a{H}_n^{t_1}$(17) where $\begin{array}{l}{C}_p^{t_0}=Q_{n-1}^{t_0}+C_a{H}_{n-1}^{t_0}-\mathrm{\Delta }t{Q}_{n-1}^{t_0}|Q_{n-1}^{t_0}|\end{array}f/2DA$, with $\begin{array}{c}{C}_a=gA/a\end{array}$.

7. Since at the same time, the flow field information at the 1D-3D coupling interface should be equally reflected on both sides, define (18) $H^{t_1}=H_n^{t_1}=H_0^{t_1}$(18) (19) $Q^{t_1}=Q_n^{t_1}=Q_0^{t_1}$(19)

8. Combine Equations (16)–(19), obtain (20) $\begin{array}{rl}{H}^{t_1}& =H_n^{t_1}=H_0^{t_1}=\frac{1}{2}\left(H_{0.5}^{t_0}-\frac{a}{gA}(Q_{0.5}^{t_0}-C_p^{t_0})\right)\end{array}$(20) (21) $\begin{array}{rl}{Q}^{t_1}& =Q_n^{t_1}=Q_0^{t_1}=C_p^{t_0}-C_a{H}^{t_1}\end{array}$(21)

9. Take the piezometric head and volumetric flow rate obtained from Equations (20) and (21) as new boundary conditions of 1D and 3D computations, and proceed to the next time step.

Figure 2. Schematic of 1D-3D coupling.

The above procedures show the case where the 1D computing domain is upstream (left) and the 3D computing domain is downstream (right). When their positions are opposite, the other characteristic equation and the Riemann invariant equation should be applied accordingly. To get good coupled results, the following points need to be noted: (1) the coupling interface should be located where the flow patterns are smooth; (2) the value of water compressibility in the 3D computation should correspond to the wave speed in the adjacent 1D computation domain; (3) at the coupling interface, the pipe’s cross-sectional area in the 1D and 3D computation should be exactly the same (note that the area of the 3D boundary may change slightly after meshing); (4) the mesh at the 3D coupling interface is preferably hexahedron with good orthogonality to ensure the volumetric flow rate of the boundary layer mesh can be calculated accurately; (5) in the case of two-phase flow, the head and flow should be corrected according to the phase fraction.

The 1D-3D coupling is implemented using the functions and classes provided by OpenFOAM, and the aforementioned 1D MOC computation programme is written in C++, and thus they can be included in the solver code of OpenFOAM and run together, allowing data to be transferred as memory variables. Within each time step loop, the code is run in the following order: (1) 3D computation code; (2) 1D computation code; (3) 1D-3D coupling code. Taking the ‘rhoPimpleFoam’ solver as an example, a brief composition of the solver code for the 1D-3D coupled algorithm is shown in Figure 3.

Figure 3. Brief composition of the solver code for 1D-3D coupled algorithm.

2.2 Additional numerical models: Gas-liquid flow and dynamic mesh

The above three basic numerical models are sufficient for general 1D-3D coupled simulations, but in practice, depending on the problem, additional numerical models may be required. For the long-distance water conveyance system in this study, gas-liquid flow phenomenon exists in the gate chamber, and the opening and closing of the gate cause dynamic boundaries; therefore, the corresponding gas-liquid flow and dynamic mesh models are also required.

2.2.1 Gas-liquid flow

The most popular multiphase flow model nowadays is VOF, where the proportion of a phase in a mesh cell is represented by a volume fraction variable, α, which takes values between 0 and 1. For gas-liquid flow, assuming that liquid is the first phase, α = 1 denotes liquid, α = 0 denotes gas, and 0 < α < 1 denotes gas-liquid mixture. The density and viscosity of the mixture are respectively expressed as (Fan & Anglart, 2020) (22) $\begin{array}{rl}\rho & =\alpha {\rho }_l+(1-\alpha ){\rho }_g\end{array}$(22) (23) $\begin{array}{rl}\mu & =\alpha {\mu }_l+(1-\alpha ){\mu }_g\end{array}$(23) where ρ, μ and α are the density, viscosity and volume fraction of liquid, respectively; the subscripts l and g denote liquid and gas, respectively. The governing equation for the volume fraction is as follows (Piscaglia et al., 2019). (24) $\frac{\mathrm{\partial }\alpha }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\alpha \mathit{U})+\mathrm{\nabla }\cdot ((1-\alpha )\alpha {\mathit{U}}_{\mathit{c}})=0$(24) where U is the flow velocity vector, U_c is the compression velocity at the intersection between the two phases; $\mathrm{\nabla }\cdot ((1-\alpha )\alpha {\mathit{U}}_{\mathit{c}})$ is the counter-gradient transport term, which guarantees conservation and boundedness.

The VOF algorithm has been implemented in the two-phase flow solvers of OpenFOAM. If both phases are incompressible, the ‘interFoam’ solver is chosen; if at least one of the phases is compressible, as in the case of water hammer, the ‘compressibleInterFoam’ solver should be chosen.

2.2.2 Dynamic mesh

In the finite volume method, when the mesh boundary is moving, the integral form of the general governing equation for a control Volume is (Demirdžić & Perić, 1990) (25) $\begin{array}{rl}& \frac{\mathrm{\partial }}{\mathrm{\partial }t}{\int }_V\rho \varphi dV+{\int }_S\rho (\mathit{U}-{\mathit{U}}_b)\varphi \cdot d\mathit{S}-{\int }_S{\mathrm{\Gamma }}_{\varphi }\mathrm{\nabla }\varphi \cdot d\mathit{S}\\ & ={\int }_V{s}_{\varphi }dV\end{array}$(25) where V is the mesh volume; S is the surface area of the mesh boundary; ϕ is the general variable; ρ is the fluid density; U is the flow velocity vector; U_b is the moving velocity vector of the mesh boundary; Γ_ϕ is the diffusion coefficient of ϕ; s_ϕ is the volumetric source of ϕ. The change rate of the mesh volume and the moving velocity of the mesh boundary satisfy the space conservation law (SCL) (26) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}\underset{V}{\int }dV-\underset{S}{\oint }{\mathit{U}}_b\cdot d\mathit{S}=0$(26) The conservation equations for mass, momentum, energy or other scalar quantities in the case of dynamic mesh can be obtained by deforming the general governing equation (Equation (25)).

Many different types of dynamic mesh have been implemented in OpenFOAM, which can be classified into mesh deformation and topology changes (Jasak, 2009). For the long-distance water conveyance system in this study, the mesh changes greatly during the opening and closing of the flat gate and radial gate; thus topology changes are needed. The most suitable dynamic mesh solver for gate motion is ‘displacementLayeredMotion’, which can be combined with layerAdditionRemoval (LAR) to make the fluid mesh under the gate increase layers during opening and decrease layers during closing, as shown in Figure 4. Note that in OpenFOAM-v2112, the above dynamic mesh solver cannot achieve the circumferential motion of the radial gate, and therefore we have modified it by using cylindrical coordinate system conversion. At the interface between the moving mesh domain and the adjacent mesh domains, the ‘cyclicACMI’ boundary is specified, enabling data exchange between the mesh domains when the size or position of the interface changes.

Figure 4. Schematic of mesh moving process for the radial gate.

3. Validation of numerical models

In this section, three water hammer cases of a reservoir-pipe-valve system are used to verify the accuracy of the 1D-3D coupled algorithm for water hammer simulations. A dynamic gate opening simulation is compared with an experiment in the literature to verify the accuracy of the gas-liquid flow simulation with dynamic mesh.

3.1 Simulations of water hammer in a reservoir-pipe-valve system

Figure 5 shows the layout of the reservoir-pipe-valve system. The total head of the upstream reservoir is H_u = 100 m; the length and diameter of the main pipeline are L = 500 m and D = 2 m, respectively; the wave speed is a = 1000 m s⁻¹; and the total head of the downstream reservoir, which is connected with the valve through a short pipe, is H_d = 0 m. The time taken for a water hammer wave to propagate one round trip is defined as one phase (2L / a). Three typical water hammer cases are simulated: the direct water hammer (the closing time is less than one phase), the first-phase water hammer (the maximum pressure at the valve occurs at the end of the first phase), and the last-phase water hammer (the maximum pressure at the valve occurs in a phase after the first phase). The simulation conditions given for each case are shown in Table 1, where the closing time is the time taken for the relative valve opening (flow area ratio) to change from 1 to 0.

Figure 5. Layout of the reservoir-pipe-valve system.

Table 1. Simulation conditions for the three typical water hammer cases.

Water hammer case	Initial discharge (m^3·s^−1)	Friction loss coefficient	Closing time (s)
Direct	20	0.02	0.5
First-phase	5	2e-5	10
Last-phase	20	0.02	10

Two algorithms are used in the simulations for comparison: (1) 1D-3D coupled algorithm, where the system is divided into two halves from the middle of the main pipe, with the upstream and downstream halves simulated using the 3D CFD method and 1D MOC, respectively; (2) pure 1D MOC. The comparison of the results is shown in Figure 6, which shows that the results of the two algorithms generally agree well, with consistent extremums, periods and trends. However, in terms of attenuation, the 1D-3D coupled algorithm is faster and closer to the experiment (S. Yang et al., 2017). This is because the 3D CFD method can capture the change in flow velocity profile along the pipe radius, which implicitly includes the influence of unsteady friction, while the steady friction model used in the 1D MOC underestimates damping (which can be improved by the unsteady friction model (Brunone et al., 2000; Urbanowicz et al., 2021)).

Figure 6. Head histories on the upstream side of the valve for the reservoir-pipe-valve system: (a) direct water hammer; (b) first-phase water hammer; (c) last-phase water hammer.

3.2 Simulation of gas-liquid flow after quickly opening a gate

According to the experimental conditions in the literature (Lobovský et al., 2014), a tank with dimensions of 1610 × 600 × 150 mm is modelled for simulation, the gate position is 600 mm from the right wall of the tank, the initial water surface height is set at 300 mm, and the gate opening speed is 3.46 m s⁻¹. The chosen OpenFOAM solver is ‘compressibleInterDyMFoam’, which is suitable for compressible two-phase flow with dynamic mesh. The type and details of the dynamic mesh have been described above. The comparison of the simulated and experimental results in the literature is shown in Figure 7, where the flow patterns are almost identical at all corresponding moments.

Figure 7. Comparison of the simulated and experimental results (Lobovský et al., 2014) for gas-liquid flow after quickly opening a gate.

4. Practical engineering application for an ultra-long water conveyance system

The 1D-3D coupled simulation method was applied to a large water conveyance project in China. The project is a pressurised gravity flow system with a total length of 194 km and a rated discharge of 265 m³ s⁻¹, as demonstrated in Figure 8. In its middle-lower reach, control gates are set to control the opening and closing of the system. The gate chamber has two layers, as shown in Figure 9, with two flat gates in the lower layer and two radial gates in the upper layer, which are bilaterally symmetrical. In front of the gate, there is an overflow surge shaft combined with the flat gate shaft, and behind the gate, there is a stilling basin and an air vent. Due to the complex structure of the gate chamber and the dynamic gas-liquid flow patterns, this gate chamber is simulated using the 3D CFD method for accuracy, while the remaining pipeline system composed of pipes, reservoirs, surge shafts and pipe connectors, is simulated using the 1D MOC, and the two domains are coupled together.

Figure 8. Schematic of the water conveyance system.

Figure 9. 3D geometry of the gate chamber fluid domain.

4.1 Computational models

4.1.1 1D computing domain (pipeline system)

The pipeline system is laid out as shown in Figure 8, where the vertical lines represent surge shafts. The total length of the pipeline system is 193.5 km, with 163.7 km upstream and 29.8 km downstream. Each single pipeline length ranges from 48 to 16273 m, and the pipe diameter ranges from 9.7 to 10.8 m. The range of the pipe’s friction loss coefficient is 8.116e−3 to 8.412e−3. The wave speed is uniformly taken as 1200 m s⁻¹. For the upstream reservoir, the normal and lowest water levels are 173.3 and 143.3 m, respectively; while for the downstream reservoir, the normal level is 86.5 m, and the highest level is  89.5 m.

4.1.2 3D computing domain (gate chamber)

Figure 9 illustrates the 3D geometry of the gate chamber, which is the fluid domain. As the gate chamber is symmetrical, to save computational resources and time, only half of the domain is simulated. The mixed mesh of tetrahedron and hexahedron is used. The upstream and downstream coupling sections are meshed with hexahedrons to ensure accurate flux for coupling, while the rest is meshed with tetrahedrons due to its complex structure, as shown in Figure 10.

Figure 10. Mesh of the 3D computing domain.

Verification of mesh independence is done. Four sets of meshes with 37, 61, 113 and 173 million cells are used to simulate two typical steady operating conditions: (1) the flat gate is half open and the radial gate is fully closed; (2) the flat gate is fully closed and the radial gate is open in small opening. The total discharge results computed by different meshes are shown in Table 2, with little difference. Considering the limited computational resources and time, the mesh with a moderate number of cells (1.13 million) is finally used for simulation, which ensures high accuracy and low cost.

Table 2. Mesh independence verification results.

	Operating condition 1		Operating condition 2
Cells number (million)	Discharge (m^3·s^−1)	Percentage (%)	Discharge (m^3·s^−1)	Percentage (%)
37	118.75	104.6	64.13	100.9
61	111.55	98.2	64.40	101.3
113	113.56	100.0	63.55	100.0
173	115.70	101.9	63.71	100.3

Operating condition 1: the flat gate is half open and the radial gate is full closed.

Operating condition 2: the flat gate is full closed and the radial gate is in small opening.

Percentage represents the ratio of each mesh to the third mesh.

4.2 Computation settings

The solver ‘compressibleInterDyMFoam’ is modified for the 1D-3D coupled simulation, and the solution of the energy equation is blocked. The water compressibility has been described in Section 2.1.2. The time step is fixed to 0.02 s to ensure stable and fast computation. The ‘RAS-realizableKE’ turbulence model is selected. The selections of discretization schemes, as shown in Table 3, are based on the tutorials in OpenFOAM and our experience in preliminary computations, balancing accuracy and stability. The temporal and momentum terms are discretized in the first-order accurate ‘Euler’ scheme and the second-order accurate ‘linearUpwind’ scheme, respectively. The laplacian scheme is ‘linear limited corrected 0.33’. In terms of turbulence, the convection terms use the ‘upwind’ scheme with first-order accuracy, and the diffusion term uses the ‘linear’ scheme with second-order accuracy. For the solution of fluid phase fraction, the two terms div(phi, alpha) and div(phirb, alpha) applies the ‘vanLeer’ scheme and ‘linear’ scheme, respectively.

Table 3. Discretization schemes (fvSchemes).

Dictionary name	Keyword	Value
ddtSchemes	default	Euler
gradSchemes	default	Gauss linear
divSchemes	div(phi,alpha)	Gauss vanLeer
	div(phirb,alpha)	Gauss linear
	div(rhoPhi,U)	Gauss linearUpwind grad(U)
	div(rhoPhi,k)	Gauss upwind
	div(rhoPhi,epsilon)	Gauss upwind
	div(((rho*nuEff)*dev2(T(grad(U)))))	Gauss linear
	div(phi,p)	Gauss linear
	div((phi + meshPhi), p)	Gauss linear
	div(phi,thermo:rho.water)	Gauss linear
	div(phi,thermo:rho.air)	Gauss linear
laplacianSchemes	default	Gauss linear limited corrected 0.33
interpolationSchemes	default	linear
snGradSchemes	default	limited corrected 0.33

The settings for solving matrix equations in the file fvSolution are as follows. The U, k and epsilon fields are solved using ‘smoothSolver’, and the p_rgh field is solved using ‘GAMG’ solver. The alpha.water field is solved using ‘smoothSolver’, with the MULESCorr option enabled. The convergence criterion for each field solution is 1e−7. For the PIMPLE algorithm settings, the momentumPredictor and correctPhi options are turned off, nOuterCorrectors is 5, nCorrectors is 3, nNonOrthogonalCorrectors is 1, and the convergence criterion for the entire flow field in one time step is 1e−3.

4.3 Boundary conditions and operating conditions

Table 4 shows the boundary types for the main boundaries of the 3D computing domain. The volumetric flow rate at the upstream inlet is a given value, specifying the corresponding velocity for U, while the water head at the downstream outlet is a given value and the corresponding pressure is specified for p_rgh. The total pressure at the boundary connected to the atmosphere is specified as 0. The gate motion is prescribed in the pointDisplacement field, with the dynamic mesh described in Section 2.2.2.

Table 4. Boundary conditions for the main boundaries.

Field name	Inlet (upstream)	Outlet (downstream)	Atmosphere
U	fixedValue	pressureInletOutletVelocity	pressureInletOutletVelocity
p_rgh	fixedFluxPressure	fixedValue	totalPressure
k	turbulentIntensityKineticEnergyInlet	inletOutlet	inletOutlet
epsilon	turbulentMixingLengthDissipationRateInlet	inletOutlet	inletOutlet
alpha.water	variableHeightFlowRate	variableHeightFlowRate	inletOutlet

The transient operating conditions simulated for this project are shown in Table 5, where each condition has two cases with fast and slow gate motions. Since there are so many conditions, only two typical ones are demonstrated in this paper: Close-1 and Open-1. When the system is to be opened, the upper radial gate is opened first to make the upstream water level drops to a certain height, and then the lower flat gate is opened; when the system is to be closed, the lower flat gate and upper radial gate are closed successively. The maximum openings of the flat gate and radial gate are 10 and 6 m, respectively. The flat gate moves at a constant speed, while for the radial gate, due to the non-linear relationship between its opening and discharge, its moving speed is divided into two stages. When the opening is below 20% (1.2 m), it moves slowly, and in the other stage (1.2–6.0 m), it moves quickly.

Table 5. The transient operating conditions simulated.

Operating condition	Description	Motion time	Upstream reservoir level	Downstream reservoir level	Key parameters of interest
Close-1	Normal closing at maximum head difference	34 min	173.3 m	86.5 m	pressure_front_max, level_surge_max, pressure_behind_min, level_stilling_min
Close-2		60 min
Open-1	Normal opening at high downstream level	34 min	173.3 m	89.5 m	pressure_behind_max, level_stilling_min
Open-2		60 min
Open-3	Normal opening at low upstream level	34 min	143.3 m	86.5 m	pressure_front_min
Open-4		60 min
Acdt-1	Emergency accident closing	1 min	173.3 m	86.5 m	level_surge_max
Acdt-2		2 min

max – maximum; min – minimum; front – in front of the gate chamber; behind – behind the gate chamber; surge – surge shaft; stilling – stilling basin.

4.4 Simulation results

4.4.1 Close-1 condition

Close-1 condition starts with the system fully open and ends with the system fully closed. The water levels of the upstream and downstream reservoirs are 173.3 and 86.5 m, respectively. The total closing time is 34 min, with 0–10 min to close the lower flat gate and 10–34 min to close the upper radial gate. A mesh with 1.13 million cells is used for the simulation, and 20 cores of the CPU Intel Xeon Gold 6248R are used for parallel running. A total of 150 min of simulation took 332 h.

4.4.1.1 Head and discharge at the inlet and outlet of the gate chamber

The histories of the water head and discharge at the inlet and outlet of the gate chamber are shown in Figure 11, in which the pure 1D MOC results are also presented for comparison. On the whole, the results of the 1D-3D coupled algorithm and the pure 1D algorithm are in good agreement, with similar trends and extremums. There is a large difference in the period of 7–20 min because the pure 1D algorithm neglects the influence of flow patterns. With the same pressure difference between the front and back of the gate, for the pure 1D algorithm, the discharge through the gate is determined only by the gate flow coefficient corresponding to the opening, whereas for the 1D-3D coupled algorithm the flow is not only related to the opening but also to the flow patterns. As can be seen from the flow patterns below, during the period of 7.5–17.5 min, the upper gate hole has an unpressurised discharge state, independent of the opening, meaning the gate flow computed by the 1D algorithm is not accurate.

Figure 11. Histories of the water head and discharge at the inlet and outlet of the gate chamber under Close-1 condition: (a) water head; (b) discharge.

4.4.1.2 Flow patterns and water levels in the gate chamber

The variation of the flow patterns and water levels in the gate chamber is shown in Figure 12. At the initial moment (0 min), all gates are open and the system is in a steady state with only the lower layer discharging. Then the flat gate starts to close. At about 7.5 min, as the lower flat gate is closed to a small opening, the water level in front of the gate rises to the bottom of the upper gate hole, and the upper layer also begins to discharge, with a flow pattern of unpressurised discharge. The lower flat gate is fully closed and the upper layer becomes the only flow passage at 10 min, when the upper radial gate begins to close. At about 17.5 min, the radial gate contacts the water flow and the water level in front of the gate starts to rise, making the flow patterns change to pressurised discharge. During 22–34 min, the upper radial gate is slowly closing with a small opening and is fully closed at 34 min. In this period, the water level in front of the gate is high and the opening is small, resulting in a high flow velocity under the gate. After 34 min, all gates are fully closed, but the water levels are still fluctuating due to the water hammer. The water level in the stilling basin reaches a minimum at about 40 min, slightly lower than the top of the downstream tunnel, but soon rebounds, and thus there is no significant intake of air into the downstream pipeline system. At about 60 min, the water in the surge shaft reaches the highest level, but no overflow. After that, the water levels still fluctuate considerably but tend to a stable state.

Figure 12. Variation of the flow patterns and water levels in the gate chamber under Close-1 condition.

4.4.1.3 Head distributions along the pipelines

The water head distributions along the entire system under the conditions of steady state and transient process conditions are shown in Figure 13, where the elevation represents the pipe top elevation. For this closing condition, the focus should be on the maximum pressure on the upstream side of the gate and the minimum pressure on the downstream side, which occur at the inlet of the gate chamber (117.22 mH2O) and the outlet of the gate chamber (−1.49 mH2O), respectively. The maximum pressure is greater than 1.3 times the maximum working pressure, and the minimum pressure is negative, which does not meet the safety requirements of the project. Countermeasures need to be taken, such as adjusting the gate motion time and law, modifying the pipeline elevation, etc.

Figure 13. Water head distributions along the system under Close-1 condition.

4.4.2 Open-1 condition

Open-1 condition starts with the system fully closed and water still, and ends with the system fully open. The water levels of the upstream and downstream reservoirs are 173.3 and 89.5 m, respectively. The total opening time is 34 min, with 0–24 min to open the upper radial gate and 24–34 min to open the lower flat gate. A mesh with 1.13 million cells is used for the simulation, which occupies 20 cores of the CPU Intel Xeon Gold 6248R for parallel running. A total of 100 min of simulation took 182 h.

4.4.2.1 Head and discharge at the inlet and outlet of the gate chamber

The histories of the water head and discharge at the inlet and outlet of the gate chamber are shown in Figure 14. The overall trends of the results of the two algorithms are consistent, with good agreement in the later stage, but there is a significant difference within the initial 30 min, mainly due to the influence of the flow patterns. In the 1D-3D coupled algorithm, the upper layer undergoes multiple flow pattern transitions. The flow pattern is initially pressurised free outflow, and then becomes pressurised submerged outflow, and finally changes to unpressurised discharge. These transitions are not considered by the pure 1D algorithm, resulting in the difference.

Figure 14. Histories of the water head and discharge at the inlet and outlet of the gate chamber under Open-1 condition: (a) water head; (b) discharge.

4.4.2.2 Flow patterns and water levels in the gate chamber

The Variation of the flow patterns and water levels in the gate chamber is shown in Figure 15. Initially, all gates are closed, and the water levels in front of and behind the gates are equivalent to those in the upstream and downstream reservoir, respectively. At 0 min, the upper radial gate begins to open. Within 0–12 min, the radial gate is in a slow opening process with a small opening, and the flow velocity under the gate is high. After 12 min, the radial gate opens faster. During 15–20 min, because of the high water level behind the gate, the flow pattern changes from the initial pressurised free outflow to pressurised submerged outflow. The water level in the stilling basin reaches its maximum at about 17.5 min, only slightly lower than the top of the stilling basin; therefore, the chamber roof may need to be raised in design. After 20 min, the bottom of the radial gate separates from the water flow due to the decrease in the water levels, resulting in the flow pattern of unpressurised discharge. At 24 min, the lower flat gate begins to open, and hence the water level in front of the gate begins to rapidly decrease. From 28 min, no water passes through the upper layer. After 34 min, all gates are fully open, and the water levels still fluctuate until 70 min.

Figure 15. Variation of the flow patterns and water levels in the gate chamber under Open-1 condition.

4.4.2.3 Head distributions along the pipelines

The water head distributions along the entire system under the conditions of steady state and transient process conditions are shown in Figure 16. For this opening condition, the maximum pressure on the downstream side of the gate should be focused, which occurs at the outlet of the gate chamber (35.14 mH2O). The maximum pressure is less than 1.3 times the maximum working pressure, meeting the safety requirements of the project.

Figure 16. Water head distributions along the system under Open-1 condition.

4.5 Flow pattern analysis of gate-controlled discharge

As can be seen from the water head and discharge at the inlet and outlet of the gate chamber mentioned above, in some cases, there are significant differences between the results of the 1D-3D coupled algorithm and the pure 1D algorithm. The essential reason is that for gas-liquid flow such as the gate-controlled discharging flow, the flow pattern is a very important factor in determining the discharge. The 1D-3D coupled algorithm, in which the local flow field is simulated by the 3D CFD method, can capture the transitions of flow patterns, whereas the pure 1D algorithm using 1D valve model to simulate gate flow neglects the influence of flow patterns.

For gate-controlled discharging flow, depending on the gate opening and the water levels in front of and behind the gate, the following various flow patterns may occur. (1) Pressurised free outflow occurs when the opening is small with a high front level and low back level, e.g. within 0–15 min under Open-1 condition (Figure 17(a)), in which the discharge is mainly determined by the opening and front level. (2) Pressurised submerged outflow occurs when the opening is small with a high back level, e.g. within 15–20 min under Open-1 condition (Figure 17(b)). In this case, the discharge is related to the opening, the back level, and the difference between the front and back levels. (3) Unpressurised discharge occurs when the opening is large and the front and back levels are both low, e.g. within 20–28 min under Open-1 condition (Figure 17(c)), where the discharge is independent of the opening. During the transient process, as the opening and the front and back level change, the above three flow patterns may occur. The determinants of discharge in each flow pattern are different and in some cases vary greatly. Therefore, when simulating gate-controlled discharging flow, it is important to consider the influence of its flow patterns.

Figure 17. Various flow patterns for gate-controlled discharging flow: (a) pressurised free outflow; (b) pressurised submerged outflow; (c) unpressurised discharge.

In practical operation, considering that the submerged discharge may induce gate vibration (especially when the submerged hydraulic jump occurs) (Erdbrink et al., 2014; Shen et al., 2018; Smok et al., 2022), it should be avoided as much as possible. For this project, under Open-1 condition, pressurised submerged outflow occurs within 15–20 min, indicating that the radial gate is opened too fast or too early. To avoid this, the slow opening stage of the radial gate (0–12 min) should be extended, so that the water levels in front of and behind the gate drop lower before quickly opening the gate.

5. Conclusions

This paper presents an efficient and accurate 1D-3D coupled simulation method based on the open source software OpenFOAM. The required numerical models and their implementation methods, including 1D MOC, water compressibility for 3D CFD, and 1D-3D coupling are described. The coupled method is applied to the simulation of the transient processes of a pressurised gravity flow water conveyance system with ultra-long pipelines and large discharge. The vivid flow patterns in the gate chamber and water head distributions along the pipelines are demonstrated. The main conclusions are as follows.

The main idea of the 1D-3D coupled algorithm in this paper is to use flow field information from the 1D and 3D computing domains in the current time step to construct their respective equations and solve them to obtain the shared boundary condition at the coupling interface for the next time step. The key and difficulty lie in constructing the equation of the 3D computing domain, which is done in this study with the Riemann invariants. Some notes on the implementation have been described. The algorithm is efficient and stable because it uses computer memory to transfer data and the calculation is explicit, which requires no iteration.

The 1D-3D coupled simulation method has the following advantages: (1) the results are more realistic and credible because the influence of flow patterns is considered; (2) the costs of computational resources and time are greatly reduced compared to the full 3D CFD method; (3) more types of flow field information can be obtained for other simulations and analyses, such as fluid-structure interaction. The method provides a means for fine simulation of the transient processes of hydraulic systems including water conveyance systems, pumped-storage power stations, etc., supporting the design and optimisation of related projects.

For gate-controlled gas-liquid flow, the influence of flow patterns on simulation results cannot be ignored. During the opening and closing process of the gate, the flow pattern may change between the following types: (1) pressurised free outflow, (2) pressurised submerged outflow, and (3) unpressurised discharge. The determinants of discharge in each flow pattern are different, and thus using the 1D valve model to simulate gate-controlled flow is inaccurate. The local 3D CFD method can capture the flow patterns and their transition, making the simulation results reliable. Therefore, for gate-controlled flow simulation, the 1D-3D coupled method is recommended.

This paper focuses on the implementation and application of the method, but the following shortcomings still need to be improved: (1) the 3D mesh is coarse and the boundary layer is not considered; (2) the unsteady friction model is not used in the 1D MOC algorithm. In this study, the coupling interface is fixed in a place with a good flow pattern and almost no gas (alpha.water > 0.9); however, in some special cases (e.g. when the gates are closed too fast), the free surface (gas-liquid interface) will move to the coupling interface, which may lead to computational errors or even divergence. Hence, movable coupling interface that can be adjusted to the position of free surface is the future research direction.

Acknowledgements

This study was sponsored by the National Natural Science Foundation of China (NSFC) (Nos. 51839008 and 51909226). The numerical simulations in this study were conducted on the supercomputing system in the Supercomputing Center of Wuhan University.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This study was sponsored by the National Natural Science Foundation of China (NSFC) (Nos. 51839008 and 51909226).