General one-dimensional hydrodynamic model flow boundary numerical value processing method

Abstract

The invention discloses a universal one-dimensional hydrodynamic model flow boundary numerical value processing method. The method includes: acquiring river basic data; carrying out finite volume element discretization on the river channel, and defining a boundary element; initializing the model; acquiring a time step dt and an outer boundary condition value at the moment t; calculating the numerical flux of the internal interface of each unit and the external interface of the water level lower boundary unit at the t moment; updating hydraulic element values of other calculation units except the flow upper boundary calculation unit; and the hydraulic element value of the flow upper boundary calculation unit is updated to the t + dt moment through a method of adding a mass source item and a momentum source item. According to the method, the one-dimensional flow boundary is processed by adding a mass source item and a momentum source item when a mass equation and a momentum equation of a flow upper boundary calculation unit are solved. No matter what flow state the inflow boundary is, only flow process data needs to be provided, and a unified numerical method is adopted for processing flow boundary conditions under all flow states, so that the universalization processing of the flow boundary is realized.

Description

Technical Field

[0001] The present invention relates to the field of water conservancy projects, in particular to the field of flood control and disaster reduction, and specifically to a general one-dimensional hydrodynamic model flow boundary numerical processing method. Background Art

[0002] One-dimensional hydrodynamic models are widely used in the field of water conservancy engineering. They can be used to calculate the real-time evolution of river floods to provide technical support for flood prevention and disaster reduction, calculate river water surface lines for water conservancy planning, and load sediment or water quality modules to expand into other application areas. Since the 1970s, with the development of computer technology and numerical analysis techniques, one-dimensional hydrodynamic model technology has continued to develop and mature. After entering the 21st century, the Godunov-type format, which originated in the field of aerodynamics and is based on solving the Riemann approximation, has the function of capturing shock waves and can simulate the movement of large gradient water surfaces. Therefore, it has gradually developed into the mainstream calculation format for one-dimensional hydrodynamic models. This type of hydrodynamic model is called a new generation of hydrodynamic models. Commonly used Godunov-type formats include the HLL format and the Roe format.

[0003] The new generation of hydrodynamic models differs significantly from the previous generation of classic implicit finite difference methods in their approach to boundary problems. These models generally rely on characteristic line theory to solve for the numerical flux at the interface of the computational unit. The original finite difference method directly solves matrices, requiring minimal processing of boundary conditions. However, using characteristic line theory to address boundary problems is significantly more complex than the original matrix-based approach, especially when dealing with flow boundary problems. This is because the inlet flow regime constantly switches between rapid, slow, and critical flow. The specific characteristic line equations used for each flow regime are different, and the basic data requirements for the flow boundary are also different. If the inlet flow regime is known to be rapid, using characteristic line theory to address the flow boundary requires not only flow process data but also corresponding water level process data. However, for slow and critical flow regimes, only flow process data is required. In practical engineering applications, flow boundary treatment presents two major challenges. First, in engineering applications, flow boundaries are typically provided through hydrological models or design flood processes, making it difficult to simultaneously provide water level data that aligns with the flow process. Second, it is difficult to determine in advance whether a river's flow boundary represents a rapid or slow flow regime; this determination can only be made during the calculation process based on the specific hydraulic element values at the current moment. In areas with steeper slopes, such as small and medium-sized rivers and mountain torrents, where rapid flow regimes are common, calculations can easily fail due to poor flow boundary treatment. In summary, although characteristic line theory is relatively mature, the treatment requirements for flow boundaries vary under different flow regimes, and no universal method has been implemented for engineering applications. Therefore, it is crucial to develop a universal numerical method for treating flow boundaries in one-dimensional hydrodynamic models. Universality means that regardless of the flow regime at the river inlet, only flow process data are required as boundary condition inputs, and the same universal numerical treatment method can be used. Summary of the Invention

[0004] The purpose of the present invention is to provide a universal flow boundary numerical processing method for the new generation of one-dimensional hydrodynamic models based on the Godunov format. Regardless of the flow state at the inlet, only the flow process data needs to be provided as the boundary condition input data, and a universal numerical processing method is used to process the flow boundary.

[0005] The present invention is achieved through the following technical solutions:

[0006] A general one-dimensional hydrodynamic model flow boundary numerical processing method is proposed. The method realizes the universal processing of one-dimensional flow boundary by adding mass source terms and momentum source terms when solving the mass equation and momentum equation of the flow upper boundary calculation unit. Regardless of whether the flow boundary is a slow flow, a rapid flow, or a critical flow state, only the flow inflow process data needs to be provided. Under rapid flow conditions, no other supporting data such as water level process is required. The specific steps of the method are as follows:

[0007] 1) Obtain one-dimensional cross-sectional data of the river channel and obtain the upstream and downstream boundary condition data required for calculation;

[0008] 2) Discrete one-dimensional river channel, the location of each section is the center of the calculation unit, the unit adjacent to the boundary condition is defined as the boundary calculation unit, the upper boundary is given as the flow boundary, without loss of representativeness, and the lower boundary is given as the water level boundary. The calculation unit corresponding to the flow boundary is the flow upper boundary calculation unit, and the calculation unit corresponding to the water level boundary is the water level lower boundary calculation unit; the outer interface of the flow upper boundary calculation unit and the outer interface of the water level lower boundary calculation unit are defined as external interfaces, and the interfaces of other calculation units are defined as internal interfaces;

[0009] 3) The hydraulic element variables (water flow area A and flow rate Q) and roughness parameters of each calculation unit are defined at the center of the calculation unit, and the initial values of the hydraulic elements and roughness parameters of each calculation unit are assigned;

[0010] 4) Get the time step dt; get the flow value Q imposed on the upper boundary calculation unit by the flow boundary at the current time t 上 The water level value Z imposed on the lower boundary calculation unit by the water level boundary 下 ;

[0011] 5) Calculate the numerical flux through the internal interface of each calculation unit and the external interface of the calculation unit through the lower boundary of the water level at time t, and the water surface gradient term and friction term of the calculation unit;

[0012] 6) Calculate the mass conservation equation and momentum conservation equation after discretization of all calculation units except the calculation unit with upper flow boundary to calculate the hydraulic element value of each calculation unit except the calculation unit with upper flow boundary at time t+dt;

[0013] 7) Add the flow boundary value Q to the discretized mass conservation equation and momentum conservation equation of the flow upper boundary calculation unit 上 The mass source term and momentum source term brought about by the equation are as follows:

[0014]

[0015] in, is the scalar form of the hydraulic element value of the upper boundary calculation unit of the flow at time t+dt, is the scalar form of the hydraulic element value of the upper boundary calculation unit of the flow at time t, dt is the time step, dx1 is the length of the upper boundary calculation unit of the flow, F 1+1 / 2 is the numerical flux of the internal interface of the flow upper boundary calculation unit, S1 is the water surface gradient term and friction term of the flow upper boundary calculation unit, is the mass source term of the calculation unit on the upper boundary of the flow, is the momentum source term of the flow upper boundary calculation unit, Q1 is the flow value of the flow upper boundary calculation unit at time t, and A1 is the water flow area of the flow upper boundary calculation unit;

[0016] Solving the above formula (1) can obtain the hydraulic element value of the flow upper boundary calculation unit at time t+dt;

[0017] 8) Let t = t + dt and repeat steps 4) to 7) until the calculation is completed.

[0018] For further optimization, in step 1), one-dimensional cross-sectional data and upstream and downstream boundary condition data of the river channel are obtained by surveying and mapping means; the distance between two adjacent cross-sectional areas is not greater than 500m, and the difference between the starting points of two adjacent measuring points in the same cross-sectional area is not greater than 10m; the upstream and downstream boundary condition data include upstream flow process data and downstream water level process data.

[0019] In step 2), the finite volume method is used to discretize the one-dimensional river channel.

[0020] Furthermore, the selection of dt in step 4) satisfies the Crown number condition, as shown in formula (2):

[0021]

[0022] Where: u is the flow velocity of the calculation unit, c is the wave velocity of the calculation unit, dx is the spatial step of the calculation unit, and dt is the time step.

[0023] Furthermore, in step 5), the complete one-dimensional Saint-Venant equations are used to describe the flow in the river. The numerical flux through the internal interface of each calculation unit is calculated using the HLL numerical format according to the hydraulic element value of each calculation unit at time t; specifically, as shown in formula (3):

[0024]

[0025] Among them, x is the spatial variable, t is the time variable, and U, F, and S are the vector expressions of the variables in the equation system, as follows:

[0026]

[0027] Where: Q is the cross-sectional flow rate, A is the cross-sectional area of the water flow, Z is the water level, g is the acceleration of gravity, and J is the resistance loss along the way. The expression is J = (n 2 Q|Q|) / (A 2 R 43 ), R is the hydraulic radius, and n is the Manning roughness coefficient.

[0028] Furthermore, in step 5), the numerical flux passing through the external interface of the water level lower boundary calculation unit is calculated by the characteristic line method.

[0029] The specific operation of step 6) is to obtain the discretized mass conservation equation and momentum conservation equation of all calculation units except the flow upper boundary calculation unit to calculate the hydraulic element value of each calculation unit except the flow upper boundary calculation unit at time t+dt; the discretized mass conservation equation and momentum conservation equation are the equations obtained by discretizing the mass conservation equation and momentum conservation equation of formula (3), as shown in formula (4):

[0030]

[0031] Where: U i t+dt and U i t are the vector forms of the hydraulic element values of the i-th calculation unit at time t+dt and time t respectively; are the vector forms of the numerical flux values of the left and right interfaces of calculation unit i, dx i is the length of the i-th computational unit; dt is the time step, S i It is the vector form of the source item value of the i-th computational unit.

[0032] The advantages and beneficial effects of the present invention are: it can effectively overcome the inherent defects of using the original characteristic line-based method to calculate the numerical flux through the left interface of the flow upper boundary calculation unit. Regardless of whether the flow boundary is a slow flow, a rapid flow or a critical flow state, it is only necessary to provide the flow inflow process data. There is also no need to first determine the flow state at the inlet and then determine the specific numerical processing method. A unified numerical method can be used to process the flow boundary conditions under all flow states, which is easy to program and implement and promote. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] Figure 1 This is a flow chart of the numerical processing method for flow boundary of a general one-dimensional hydrodynamic model in Example 1 of the present invention;

[0034] Figure 2 This is a schematic diagram of a finite volume element according to embodiment 1 of the present invention;

[0035] Figure 31 is a comparison chart of the calculation results of the test example with analytical solution of Example 1 of the present invention, (a) slow flow flow boundary result, (b) rapid flow flow boundary result. DETAILED DESCRIPTION

[0036] The following is combined with Figure 1 The present invention is further described with reference to the accompanying drawings and examples.

[0037] The present invention provides a general one-dimensional hydrodynamic model flow boundary numerical processing method. The method is based on the idea of conservation of momentum and mass. By adding mass source terms and momentum source terms when solving the mass equation and momentum equation of the flow upper boundary calculation unit, the one-dimensional flow boundary is processed. This method can effectively overcome the problem of using the original characteristic line method to calculate the flow through the outer side of the flow upper boundary calculation unit (attached). Figure 2 The inherent defects of numerical flux calculation of the interface (left side in the figure) are eliminated by the method provided by the present invention. Regardless of whether the flow boundary is a slow flow, a rapid flow, or a critical flow state, only the flow inflow process data needs to be provided. Under rapid flow conditions, the supporting water level process data is no longer required. Similarly, there is no need to first determine the flow state at the inlet and then determine the specific numerical processing method. A unified numerical method can be used to process the flow boundary conditions under all flow states, realizing universal processing of the flow boundary. The method includes the following specific steps:

[0038] 1) Obtain basic river data:

[0039] One-dimensional cross-sectional data of the river channel is obtained through surveying and mapping. The distance between adjacent sections is no more than 500 meters, and the starting points of adjacent measurement points on the same section should not differ by more than 10 meters. High-efficiency surveying methods such as drones equipped with LiDAR sensors can be used to obtain laser point cloud data above the river surface. Survey vessels equipped with ADCP equipment can be used to map the subsurface topography of the river channel. The above- and below-surface topographic data are then fused, and interpolation is used to obtain specific cross-sectional data. The upstream and downstream boundary condition data required for calculations generally include upstream flow process data and downstream water level process data.

[0040] 2) Discretize the river channel into finite volume units and define the boundary units:

[0041] The finite volume method is used to discretize the one-dimensional river channel. The location of each section is the center of the calculation unit. The unit adjacent to the boundary is defined as the boundary calculation unit. The upper boundary is usually given as the flow boundary, which is representative. The lower boundary is given as the water level boundary. The calculation unit adjacent to the flow boundary is the flow upper boundary calculation unit, and the calculation unit adjacent to the water level boundary is the water level lower boundary calculation unit. The left interface of the flow upper boundary calculation unit and the right interface of the water level lower boundary calculation unit are defined as external interfaces, and the interfaces of other calculation units are defined as internal interfaces. Figure 2 shown.

[0042] 3) Model initialization:

[0043] The hydraulic element variables (water flow area A and flow rate Q) and roughness parameters of each calculation unit are defined at the center of the calculation unit, and the initial values of the hydraulic elements and the roughness parameters are assigned to each calculation unit.

[0044] 4) Get the time step dt and outer boundary condition value at the current time t:

[0045] The calculation is solved using the explicit format to obtain the calculation step length dt that meets the Crown number condition; the specific formula is as follows:

[0046]

[0047] Where: N cfl is the Crown number, u is the flow velocity of the calculation unit (m / s), c is the wave velocity of the calculation unit (m / s), dx is the spatial step size of the calculation unit (m), and dt is the time step size (s).

[0048] According to the upper boundary flow process data and lower boundary water level process data provided, the flow value Q imposed on the upper boundary calculation unit by the flow boundary at the current time t is obtained by linear interpolation. 上 The water level value Z imposed on the lower boundary calculation unit by the water level boundary 下 .

[0049] 5) Calculate the internal interface of each unit and the right side of the boundary unit below the water level at time t ( Figure 2 ) interface: The complete one-dimensional Saint-Venant equations are used to describe the flow in the river, as shown in Equation (3):

[0050]

[0051] Among them, x is the spatial variable, t is the time variable, and U, F, and S are the vector expressions of the variables in the equation system, as follows:

[0052]

[0053] Where: Q is the cross-sectional flow rate (m 3 / s), A is the cross-sectional area of water (m 2 ), Z is the water level (m), g is the acceleration due to gravity (m / s 2 ), J is the resistance loss along the way, and its expression is J=(n 2 Q|Q|) / (A 2 R 43 ), R is the hydraulic radius (m), n is the Manning roughness coefficient (m 1 / 3 / s).

[0054] Based on the hydraulic element values of each computational cell at time t, the numerical flux through the internal interface of each computational cell is calculated using the HLL numerical format. The specific HLL numerical format can be found in the following literature (Zhang Dawei, Numerical Simulation of Dam Breach Flow Based on Godunov Format [M]. China Water Resources and Hydropower Press, Beijing, 2014, 12). The numerical flux through the outer interface of the computational cell at the lower boundary of the water level is calculated using the characteristic line method. The specific theory of the characteristic line method can be found in the following literature (Toro FE, Riemann solvers and numerical methods for fluid dynamics [M]. Berlin: Springer-Verlag, 1999.).

[0055] Calculate the water surface gradient term and friction term for each calculation unit:

[0056] Calculate the water surface gradient term and friction term of each calculation unit. For the specific calculation process, please refer to the following literature (Zhang Dawei, Cheng Xiaotao, Huang Jinchi, etc. A highly applicable mathematical model for complex open channel flow [J]. Journal of Hydraulic Engineering, 2010, 41(4): 531-536.).

[0057] 6) Update the hydraulic element values of all calculation units except the upper boundary calculation unit of flow to t+dt:

[0058] The mass conservation equation and momentum conservation equation after discretization of all calculation units except the flow upper boundary calculation unit are solved to calculate the hydraulic element values of each calculation unit except the flow upper boundary calculation unit at time t+dt.

[0059] The mass conservation equation and momentum conservation equation after the discretization are the mass conservation equation and momentum conservation equation of formula (3) after the discretization. Formula (3) is numerically discretized by the finite volume method based on the central format to obtain formula (4):

[0060]

[0061] Where: U i t+dt and U i t are the vector forms of the hydraulic element values of the i-th calculation unit at time t+dt and time t respectively; are the vector forms of the numerical flux values of the left and right interfaces of calculation unit i, dx i is the length of the i-th computational unit; dt is the time step, S i It is the vector form of the source item value of the i-th computational unit.

[0062] 7) Update the hydraulic element value of the flow upper boundary calculation unit to time t+dt by adding mass source terms and momentum source terms:

[0063] Add the flow boundary value Q to the mass conservation equation and momentum conservation equation after discretization in the flow upper boundary calculation unit 上 The mass source term and momentum source term brought about by the equation are as follows:

[0064]

[0065] in, is the scalar form of the hydraulic element value of the upper boundary calculation unit of the flow at time t+dt, is the scalar form of the hydraulic element value of the upper boundary calculation unit of the flow at time t, dt is the time step, dx1 is the length of the upper boundary calculation unit of the flow, F 1+1 / 2 is the numerical flux on the right interface of the flow upper boundary calculation unit, S1 is the water surface gradient term and friction term of the flow upper boundary calculation unit, is the mass source term of the calculation unit on the upper boundary of the flow, is the momentum source term of the flow upper boundary calculation unit, Q1 is the flow value of the flow upper boundary calculation unit at time t, and A1 is the water flow area of the flow upper boundary calculation unit.

[0066] By solving the above formula (1), we can obtain the hydraulic element value of the flow upper boundary calculation unit at time t+dt.

[0067] 8) Let t = t + dt and repeat steps 4) to 7) until the calculation is completed.

[0068] Figure 3 The results are for a steady-state flow example in an open channel with a throat. The calculation area for this example is a 3-meter-long rectangular channel. In the middle of the channel, the river width shrinks, forming a throat-like channel. At the throat, the riverbed elevation also rises accordingly.

[0069] The expression for the change of river width b(x) is:

[0070]

[0071] Bottom elevation z b The expression of (x) change is as follows:

[0072]

[0073] Under different inflow flow rates, different flow patterns of slow flow and rapid flow will appear at the inlet. The specific working conditions are as follows:

[0074] (1) Slow flow rate boundary: inlet flow rate Q = 1.565m 3 / s, downstream water level Z = 1.0m.

[0075] (2) Rapid flow boundary: inlet flow Q = 5.322m 3 / s, downstream water level Z = 1.0m.

[0076] The influence of friction is not considered in the calculation. The initial water level of all calculation cells is given as 1.0m, and the initial flow rate is given as 0m. The calculation area is discretized using uniform cells with a spatial step size of 0.02m, and is divided into 150 calculation cells in total.

[0077] The flow boundary processing method proposed in the present invention is used to calculate the above two working conditions respectively, and the results are compared with the theoretical analytical solution of the example. The results of the analytical solution of the example can be referred to the following literature (Hubbard M E. On the accuracy of one dimensional models of steady converging / diverging open channel flows [J]. International journal for numerical methods in fluids, 2001, 35 (7): 785-808.). Figure 3 The comparison results of the numerical solution and analytical solution of the water surface calculation when the calculation reaches a stable state are shown in Figure 2. Figure 3 (a) is the water surface curve result under the slow flow boundary, Figure 3 (b) is the water surface curve result under the rapid flow boundary. It can be seen from the two figures that the numerical value of the water level is in good agreement with the analytical solution under the slow inflow boundary and the rapid inflow boundary, indicating that the boundary processing method proposed in this invention is successful in terms of versatility.

[0078] The above embodiments are only partial embodiments of the present invention and cannot cover the entire present invention. Based on the above embodiments and drawings, those skilled in the art can obtain more implementation methods without paying any creative work. Therefore, these implementation methods obtained without paying any creative work should be included in the scope of protection of the present invention.

Claims

1. A general one-dimensional hydrodynamic model flow boundary numerical processing method, characterized by: The following steps are involved: 1) Obtain one-dimensional cross-sectional data of the river channel and upstream and downstream boundary condition data; 2) Discrete one-dimensional river channel, the location of each section is the center of the calculation unit, the unit adjacent to the boundary is defined as the boundary calculation unit, the upper boundary is given as the flow boundary, the lower boundary is given as the water level boundary, the calculation unit adjacent to the flow boundary is the flow upper boundary calculation unit, and the calculation unit adjacent to the water level boundary is the water level lower boundary calculation unit; the outer interface of the flow upper boundary calculation unit and the outer interface of the water level lower boundary calculation unit are defined as the external interface, and the interfaces of other calculation units are defined as the internal interface; 3) Set initial values and roughness parameters for hydraulic elements of each calculation unit; 4) Obtain the time step dt and the flow value Q imposed by the flow boundary on the flow upper boundary calculation unit at the current time t 上 The water level value Z imposed on the water level lower boundary calculation unit by the water level boundary 下 ; 5) Calculate the numerical flux through the internal interface of each calculation unit and the external interface of the calculation unit through the lower boundary of the water level at time t, and the water surface gradient term and friction term of the calculation unit; 6) Calculate the hydraulic element values of each calculation unit except the calculation unit with upper flow boundary at time t+dt; 7) Solve equation (1) to obtain the hydraulic element value of the flow upper boundary calculation unit at time t+dt: in, and are the vector forms of the hydraulic element values of the upper boundary calculation unit of the flow at time t+dt and t, dx1 is the length of the upper boundary calculation unit of the flow, F 1+1 / 2 is the numerical flux of the internal interface of the flow upper boundary calculation unit, S1 is the water surface gradient term and friction term of the flow upper boundary calculation unit, is the mass source term of the calculation unit on the upper boundary of the flow, is the momentum source term of the flow upper boundary calculation unit, Q1 is the flow value of the flow upper boundary calculation unit at time t, and A1 is the water flow area of the flow upper boundary calculation unit; 8) Let t = t + dt and repeat steps 4) to 7) until the calculation is completed.

2. A general one-dimensional hydrodynamic model flow boundary numerical processing method according to claim 1, characterized in that: In step 1), one-dimensional cross-sectional data and upstream and downstream boundary condition data of the river channel are obtained by surveying and mapping means; the distance between two adjacent cross-sections is not greater than 500m, and the difference between the starting points of two adjacent measuring points in the same section is not greater than 10m; the upstream and downstream boundary condition data include upstream flow process data and downstream water level process data.

3. The method for numerically processing flow boundaries of a universal one-dimensional hydrodynamic model according to claim 1 is characterized by: Step 2) Use the finite volume method to discretize the one-dimensional river channel.

4. A general one-dimensional hydrodynamic model flow boundary numerical processing method according to claim 1, characterized in that: In step 4), the time step dt at the current time t is obtained, and the calculation is solved using the explicit format. The selection of dt satisfies the Crown number condition, as shown in formula (2): Where: N cfl is the Crown number, u is the flow velocity of the calculation unit, c is the wave velocity of the calculation unit, dx is the spatial step of the calculation unit, and dt is the time step.

5. The method for numerically processing flow boundaries of a universal one-dimensional hydrodynamic model according to claim 1, characterized in that: In step 5), the complete one-dimensional Saint-Venant equations are used to describe the flow in the river. The HLL numerical format is used to calculate the numerical flux through the internal interface of each calculation unit according to the hydraulic element value of each calculation unit at time t. The one-dimensional Saint-Venant equations are specifically shown in formula (3): Among them, x is the spatial variable, t is the time variable, and U, F, and S are the vector expressions of the variables in the equation system, as follows: Where: Q is the cross-sectional flow rate, A is the cross-sectional area of the water flow, Z is the water level, g is the acceleration of gravity, and J is the resistance loss along the way. The expression is J = (n 2 Q|Q|) / (A 2 R 4 / 3 ), R is the hydraulic radius, and n is the Manning roughness coefficient.

6. A general one-dimensional hydrodynamic model flow boundary numerical processing method according to claim 1, characterized in that: In step 5), the numerical flux through the external interface of the calculation unit at the lower boundary of the water level is calculated by the characteristic line method.

7. The method for numerically processing flow boundaries of a universal one-dimensional hydrodynamic model according to claim 5, characterized in that: The specific operation of step 6) is to obtain the discretized mass conservation equation and momentum conservation equation of all calculation units except the flow upper boundary calculation unit to calculate the hydraulic element value of each calculation unit except the flow upper boundary calculation unit at time t+dt; the discretized mass conservation equation and momentum conservation equation are the equations obtained by discretizing the mass conservation equation and momentum conservation equation of formula (3), as shown in formula (4): Where: U i t+dt and U i t are the vector forms of the hydraulic element values of the i-th calculation unit at time t+dt and time t respectively; are the vector forms of the numerical flux values of the left and right interfaces of calculation unit i, dx i is the length of the i-th computational unit; dt is the time step, S i It is the vector form of the source item value of the i-th computational unit.

Images