River flow calculation method

Abstract

The invention relates to the technical field of hydrographic survey, in particular to a riverway water flow calculation method, which comprises the following steps: acquiring topographic data, hydrological data and roughness data, and constructing a one-dimensional riverway hydrodynamic numerical model; the one-dimensional river channel hydrodynamic numerical model is based on a discontinuous finite element method and a Runge-Kutta time discrete shallow water equation set, and a numerical flux is solved by combining an essential oscillation-free ENO slope limiter and an HLL Riemann solver; aiming at the confluence situation of tributaries in different intervals, the interval flow is calculated by utilizing an iterative method, including initially assuming that the tributary flow is zero, simulating non-tributary outflow by using a one-dimensional river hydrodynamic numerical model, determining the difference between the simulated non-tributary outflow and actual site flow measurement, and adjusting the tributary flow; and iteration is repeated until the error between the simulated outlet flow and the actual outlet flow reaches the preset precision. By constructing a one-dimensional river water power numerical model based on a Runge-Kutta time discrete shallow water equation set, high-precision simulation of river water flow calculation under the condition of branch flow convergence in different intervals is realized.

Description

Technical Field

[0001] This application relates to the field of hydrological measurement technology, and more specifically, to a method for calculating river flow. Background Technology

[0002] River flow calculations are a crucial basis for reservoir flood control scheduling. Excessive calculation errors in tributaries within a given area can cause the calculated river flow to deviate from reality, severely reducing the accuracy of the calculations. Therefore, accurately determining the inflow process of tributaries within a given area is of great significance for refined flood control.

[0003] In practical engineering, the flow rate of a given interval is usually estimated proportionally based on the catchment area of ​​the watershed, which introduces significant errors. Furthermore, some studies use hydrological models or directly calculate the total inflow of the interval from the difference in water volume between the inflow and outflow of the river channel, but this is insufficient to meet the needs of real-time flood management. Other studies calculate the total inflow of the interval based on hydrodynamic models and data assimilation methods, but this method cannot derive the individual inflow process of a specific tributary. Additionally, some scholars have estimated the inflow of tributaries in the interval by directly calculating the flow from downstream to upstream, including Muskingen inversion and inverse solution of the Saint-Venant equations; however, the stability of the inversion results is poor, and the estimated inflow process has considerable uncertainty. Summary of the Invention

[0004] In view of this, this application provides a method for calculating river flow to achieve high-precision simulation of river flow under different tributary confluence scenarios.

[0005] The technical solution provided in this application is as follows:

[0006] A method for calculating river flow includes:

[0007] Topographic data, hydrological data, and roughness data are acquired to construct a one-dimensional river channel hydrodynamic numerical model. The one-dimensional river channel hydrodynamic numerical model is based on the discontinuous finite element method and the Runge-Kutta time-discrete shallow water equations, and combines the inherently non-oscillating ENO slope limiter and the HLL Riemann solver to solve for numerical flux.

[0008] For different tributary confluence scenarios, an iterative method is used to calculate the interval flow. The iterative method includes initially assuming that the tributary flow is zero, simulating the outflow without tributaries using the one-dimensional river hydrodynamic numerical model, and determining the difference between the simulated outflow without tributaries and the actual flow measured at the station. The tributary flow is adjusted according to the difference and the river hydraulic characteristics and water propagation velocity, and the iteration is repeated until the error between the simulated outflow and the actual outflow reaches the preset accuracy.

[0009] In one possible implementation, the governing equations of the model are the Saint-Venant equations for one-dimensional shallow water flow, including continuity equations and motion equations, used to describe the continuity and momentum conservation of river flow.

[0010] The continuity equation is expressed as:

[0011]

[0012] In the formula, A is the cross-sectional area of ​​the water passage; Q is the cross-sectional flow rate;

[0013] The equation of motion is expressed as:

[0014]

[0015] In the formula, S x S represents the riverbed slope. f I1 is the friction slope; I2 is the hydrostatic pressure; g is the wall pressure of the natural river channel; g is the acceleration due to gravity.

[0016] Among them, the riverbed slope S x The calculation formula is:

[0017]

[0018] In the formula, h b denoted as the riverbed elevation, and x as the distance along the river channel.

[0019] Friction slope S f This represents the energy loss caused by friction between the riverbed and the water body, calculated using Manning's formula, and is expressed as:

[0020]

[0021] In the formula, n is the Manning roughness coefficient; It is the hydraulic radius, where P is the wetted perimeter;

[0022] The formula for calculating hydrostatic pressure I1 is:

[0023]

[0024] The formula for calculating the wall pressure I2 of a natural river channel is:

[0025]

[0026] Where h is the water depth, which varies with time and space; η is the water level; B(x,η) is the cross-sectional width of the wetted perimeter of the open channel when the water level of the natural river is η.

[0027] The partial differential equation in its conserved form can be expressed as a governing equation as follows:

[0028]

[0029] In the formula, U is a vector of conserved variables. F is the flux vector. S is the source term vector.

[0030] The Jacobian matrix of the governing equations is expressed as:

[0031]

[0032] In the formula, c is the gravitational wave velocity; u is the average flow velocity; and T is the width of the water surface.

[0033] Its eigenvalues ​​and eigenvectors are represented as follows:

[0034]

[0035] In one possible implementation, after constructing the one-dimensional river hydrodynamic numerical model, the method further includes performing forward and reverse verification on the one-dimensional river hydrodynamic numerical model.

[0036] The positive verification involves comparing the simulated water depth of the one-dimensional river hydrodynamic numerical model with the theoretical water depth, and calculating the mean absolute error (MAE) and mean relative error (MRE) to evaluate the accuracy of the model.

[0037] The reverse verification is performed by calculating the water flow from downstream to upstream, comparing the simulated water depth of the one-dimensional river hydrodynamic numerical model with the theoretical water depth, and calculating the mean absolute error (MAE) and mean relative error (MRE) to evaluate the accuracy of the model.

[0038] In one possible implementation, the one-dimensional channel hydrodynamic numerical model is based on the discontinuous finite element method and the Runge-Kutta time-discrete shallow water equations, including:

[0039] The one-dimensional problem domain is divided into Ne elements and Ne+1 nodes using the discontinuous finite method. Variables within an element are continuous, while those at element boundaries are discontinuous.

[0040] The integral formula for the equation applied to each element in the discontinuous Galerkin method is as follows:

[0041]

[0042] In the formula, and These are the coordinates of the starting and ending points of the finite element element, N. i and N j U is a shape function. jLet F represent the conserved variables of flow rate Q and cross-sectional area A, F represent the flux vector, and S represent the source term vector. Represents numerical flux; N is a shape function;

[0043] The mass and momentum equations of the shallow water equations are both converted to discontinuous Galerkin form, expressed as:

[0044]

[0045] Where z represents the water level;

[0046] The discontinuous finite element method for discretizing the one-dimensional shallow water equations is expressed as:

[0047]

[0048] In the formula, Δx represents the length of the element, and ξ represents the local coordinates, which maps the physical element to ξ∈[-1,1], facilitating numerical integration using the Gaussian integral formula. N i (ξ) represents a shape function; S i (ξ) represents the source term;

[0049] In the formula, p is the mass matrix, expressed as:

[0050]

[0051] In the formula, Ω e Representing unit e, This represents the integral over the e-th unit;

[0052] Discretization in time is performed using the Runge-Kutta method, which reduces the total variation by third order, and is expressed as follows:

[0053]

[0054] In the formula, Δt represents the time step, Δt i The value represents the time step of the i-th unit. CFL stands for Courant-Friedrichs-Lewy condition, which is used to ensure the computational stability of the numerical solution. n represents the n-th step in the time-progression calculation, and the superscript indicates that the value is taken at the n-th time step.

[0055] One possible implementation combines an inherently non-oscillating ENO slope limiter and an HLL Riemann solver to solve for the numerical flux, including:

[0056] By reconstructing variables with finite slopes using adjacent elements, the equations of this element and its preceding and following adjacent elements are solved simultaneously to obtain:

[0057]

[0058] By solving the system of equations simultaneously for a and b, c and d, and determining the slope and intercept of the first-order polynomial with the smaller slope, e(i) represents the i-th unit. This represents the average cross-sectional area of ​​the water passage in the i-th unit. This represents the average flow rate in the i-th unit;

[0059] In the Galerkin method for discontinuous hydrodynamics in one-dimensional channels, the numerical flux at the boundary between elements is... The HLL flux function is used instead, and is expressed as:

[0060]

[0061] Among them, F L F represents the flux of the cell on the left side of the interface. R This is represented as the flux of the right-hand cell;

[0062] In the formula, S L The wave velocity on the left side at the boundary between elements is expressed as:

[0063]

[0064] S R The wave velocity on the right side at the element-to-element boundary is expressed as:

[0065]

[0066] Where u* and c* are respectively represented as:

[0067]

[0068] In the formula, L and R represent the limit values ​​on both sides of the element boundary, respectively; u - u+ represents the flow velocity to the left of the element boundary, and u+ represents the flow velocity to the right of the element boundary; h + h represents the water depth on the right side of the cell. - The value represents the water depth on the left side of the cell; u* represents the intermediate computational cost with the same flow velocity at the boundary; and c* represents the intermediate computational cost related to the wave velocity at the boundary.

[0069] In one possible implementation, the calculation of the interval flow includes:

[0070] The river channel is divided into multiple calculation intervals based on the location where the tributaries flow into it;

[0071] Within each calculation interval, the flow changes before and after the tributary's inflow are simulated using the one-dimensional hydrodynamic model.

[0072] For each calculation interval, the spatiotemporal distribution characteristics of the interval flow are calculated based on the simulation results;

[0073] By combining the flow changes in different calculation intervals, the relationship between interval flow and reservoir discharge flow is evaluated;

[0074] A disturbance is applied at the confluence of the tributaries, and the tributary lag time is determined by monitoring the response time of the cross section. The input time of the tributary flow is then adjusted based on the lag time to simulate the impact of the tributaries on the river flow.

[0075] In one possible implementation, the calculation of the interval flow rate further includes, depending on the location where the tributary flows in:

[0076] Calculate the inter-regional flow characteristics at the upstream, midstream and downstream confluence points, including peak flow, peak occurrence time and average flow.

[0077] By comparing the characteristics of inter-regional flow under different scenarios, the influence of the location of tributary confluence on inter-regional flow is analyzed.

[0078] One possible implementation involves iterating repeatedly until the error between the simulated and actual outlet flow rates reaches a preset accuracy, including:

[0079] The tributary flow is gradually corrected by iterative method so that the error between the simulated outlet flow and the actual observed flow is within the preset accuracy range.

[0080] In each iteration, a correction strategy is adopted to adjust the tributary flow based on the tributary inflow location, in order to adapt to the water flow changes under different scenarios and dynamically optimize the simulation results.

[0081] Compared with the prior art, the technical solution of this application has the following beneficial effects:

[0082] The technical solution of this application constructs a one-dimensional river hydrodynamic numerical model based on the discontinuous finite element method and the Runge-Kutta time-discrete shallow water equations, and combines the inherently non-oscillating ENO slope limiter with the HLL Riemann solver to achieve high-precision simulation of river flow under different tributary confluence scenarios. This scheme can effectively handle the influence of tributary confluence location on flow calculation, significantly improving the accuracy and reliability of simulation results. Dynamic correction of tributary discharge through iterative methods further optimizes simulation accuracy, ensuring a high degree of fit between simulated outflow and actual observed discharge. Furthermore, the dual verification mechanism of forward and reverse validation enhances the model's stability and generalization ability, enabling it to adapt to various complex flow conditions and tributary confluence scenarios. This effectively solves the simulation error problem caused by different tributary confluence locations in existing technologies, providing a more accurate and efficient technical means for river flow calculation. Attached Figure Description

[0083] Figure 1This is a flowchart of a river flow calculation method provided in Embodiment 1 of this application.

[0084] Figure 2 This is a flowchart of a river flow calculation method based on RGDK provided in Embodiment 2 of this application.

[0085] Figure 3 This is a schematic diagram of a one-dimensional region of a linear element in the discontinuous finite element method provided in Embodiment 2 of this application.

[0086] Figure 4 This is a comparison diagram of the initial model and theoretical water depth provided for Embodiment 2 of this application.

[0087] Figure 5 The simulation-theoretical water depth comparison diagram provided for the parameter calibration of Embodiment 2 of this application.

[0088] Figure 6 This is a schematic diagram of the simulated water depth-flow curve for a river channel provided in Embodiment 2 of this application.

[0089] Figure 7 This is a comparison chart of reverse simulation and theoretical water depth provided in Embodiment 2 of this application.

[0090] Figure 8 This is a schematic diagram of the actual inlet and outlet flow and the outflow flow without tributaries provided for Embodiment 2 of this application.

[0091] Figure 9 This is a schematic diagram of the tributary confluence section and observation section provided in Embodiment 2 of this application.

[0092] Figure 10 This is a schematic diagram of the response process of the inlet flow at the monitoring section provided in Embodiment 2 of this application.

[0093] Figure 11 This is a schematic diagram of the response process of tributary flow at the monitoring section under the confluence scenario 5km upstream provided in Embodiment 2 of this application.

[0094] Figure 12 This is a schematic diagram of the first iteration of tributary flow under the confluence scenario at 5km upstream, provided in Embodiment 2 of this application.

[0095] Figure 13 This is a comparison diagram of simulated outflow and actual outflow under the confluence scenario 5km upstream provided in Embodiment 2 of this application.

[0096] Figure 14 This is a schematic diagram of the upstream reservoir discharge and tributary flow process under the confluence scenario at 5km upstream, provided in Embodiment 2 of this application.

[0097] Figure 15This is a schematic diagram of the response process of tributary flow at the monitoring section under the confluence scenario at 50km in the middle reaches, as provided in Embodiment 2 of this application.

[0098] Figure 16 This is the first iteration of tributary flow under the confluence scenario at 50km midstream provided in Embodiment 2 of this application.

[0099] Figure 17 This is a comparison diagram of simulated outflow and actual outflow under the confluence scenario at 50km midstream provided in Embodiment 2 of this application.

[0100] Figure 18 This is a diagram illustrating the upstream reservoir discharge and tributary flow process under a confluence scenario at 50km downstream, as provided in Embodiment 2 of this application.

[0101] Figure 19 This is a schematic diagram of the response process of tributary flow at the monitoring section under the confluence scenario at a downstream location of 95km, as provided in Embodiment 2 of this application.

[0102] Figure 20 This is the first iteration of the tributary flow under the confluence scenario at a downstream location of 95km provided in Embodiment 2 of this application.

[0103] Figure 21 This is a comparison diagram of simulated outflow and actual outflow under the confluence scenario at a downstream location of 95km provided in Embodiment 2 of this application.

[0104] Figure 22 This is a schematic diagram of the upstream reservoir discharge and tributary flow process under the confluence scenario at a downstream location of 95km, as provided in Embodiment 2 of this application.

[0105] Figure 23 This is a schematic diagram of the interval flow process under three scenarios provided in Embodiment 2 of this application.

[0106] Figure 24 These are interval flow characteristic diagrams for three scenarios provided in Embodiment 2 of this application. Detailed Implementation

[0107] The technical solutions in the embodiments of this application will be clearly and completely described below with reference to the embodiments of this application. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of this application without creative effort are within the scope of protection of this application.

[0108] To address the problems mentioned in the background, this application proposes to construct a one-dimensional river hydrodynamic numerical model. The model calculates the tributary backflow flow or the confluence flow along the river into concentrated inflow flows at different locations based on the flow process observed at the river flood control section and the reservoir discharge process, thereby providing technical support for refined flood control.

[0109] Example 1

[0110] See Figure 1 This is a flowchart of a river flow calculation method provided in Embodiment 1 of this application. Figure 1 As shown, the specific implementation steps of the above method include:

[0111] Step 101: Obtain topographic data, hydrological data, and roughness data to construct a one-dimensional river channel hydrodynamic numerical model.

[0112] In this application, the one-dimensional channel hydrodynamic numerical model is a one-dimensional shallow-water equation hydrodynamic numerical model based on the RKDG scheme with an ENO-type slope limiter, which is a numerical approximation method for solving partial differential equations. The Runge-Kutta Discontinuous Galerkin method (RKDG) numerical model uses the ENO (Essentially Non-Oscillatory) type essentially non-oscillatory slope limiter technique to solve the non-physical pseudo-oscillation problem that occurs in the discontinuous finite element method. In the RKDG method, the third-order Runge-Kutta time integral with decreasing total variation is combined with the discontinuous finite element spatial discretization to improve numerical stability and accuracy.

[0113] Specifically, the one-dimensional river hydrodynamic numerical model described in this application embodiment is based on the discontinuous finite element method and the Runge-Kutta time-discrete shallow water equations, combined with the inherently non-oscillating ENO slope limiter and the HLL Riemann solver to solve for numerical flux.

[0114] Step 102: For different tributary confluence scenarios, calculate the interval flow using an iterative method. This iterative method includes initially assuming zero tributary flow, simulating tributary-free outflow using a one-dimensional river hydrodynamic numerical model, and determining the difference between the simulated tributary-free outflow and the actual flow measured at the station. Based on this difference, the tributary flow is adjusted according to the river's hydraulic characteristics and flow propagation velocity. This iteration is repeated until the error between the simulated and actual outflow reaches a preset accuracy.

[0115] Example 2

[0116] Embodiment 2 of this application is a further explanation of the technical solution provided in Embodiment 1 of this application. See also... Figure 2This is a flowchart of a river flow calculation method based on the RKDG model provided in Embodiment 2 of this application. Figure 2 As shown, the specific implementation steps of the above method include:

[0117] Step 201: Obtain topographic data, hydrological data, and roughness data.

[0118] As a feasible approach, a corresponding terrain file is created based on pre-determined river channel conditions. In this process, the river channel is divided into m cross-sections, and n nodes are set on each cross-section. Based on the division of the river channel terrain, three key terrain files are generated: S.txt, X.txt, and Y.txt. The S.txt file contains m numbers, sequentially providing the mileage of each cross-section. The X.txt file has n rows and m columns, providing the x-coordinates of the nodes on each cross-section. The Y.txt file has n rows and m columns, providing the y-coordinates of the nodes on each cross-section.

[0119] The aforementioned hydrological data includes the initial flow and water level of each cross-section, the inlet flow process and tributary flow process at each time step, which are required for the one-dimensional river hydrodynamic numerical model. The outlet is assumed to be free outflow. The corresponding hydrological files include Q_In.txt, InitialDischarge.txt, InitialWaterLevel.txt, and TributaryData_Q.txt. The Q_In.txt file provides the inlet flow process in cubic meters per second (m / s), with one data point every hour. The InitialDischarge.txt file provides the initial flow rate, containing m numbers, sequentially providing the initial flow rate at each cross-section in cubic meters per second. The InitialWaterLevel.txt file provides the initial water level, containing m numbers, sequentially providing the initial water level at each cross-section. The first line of the TributaryData_Q.txt file is the cross-section number of the tributary's location; from the second line onwards, each line shows the flow rate of each tributary at different times. Positive flow values ​​indicate inflow, and negative values ​​indicate outflow.

[0120] The roughness data above is corresponding to the roughness file ManningRoughnessCoefficient.txt, which contains m numbers, giving the Manning roughness coefficient for each cross section.

[0121] Step 202: Construct a one-dimensional river channel hydrodynamic numerical model.

[0122] Specifically, the governing equations of the aforementioned one-dimensional river hydrodynamic numerical model are the Saint-Venant equations for one-dimensional shallow water flow, including a continuity equation with the law of conservation of mass and a motion equation with the law of conservation of momentum. The continuity equation with the law of conservation of mass is expressed as:

[0123]

[0124] In equation (1), A is the cross-sectional area of ​​the water passage (m²). 2 Q is the cross-sectional flow rate (m³ / s). 3 / s).

[0125] The equation of motion that obeys the law of conservation of momentum is expressed as:

[0126]

[0127] In equation (2), S x S represents the riverbed slope. f I1 is the friction slope; I2 is the hydrostatic pressure; I3 is the wall pressure of the natural river channel; g is the acceleration due to gravity (m / s²). 2 ).

[0128] Among them, the riverbed slope S x The calculation formula is:

[0129]

[0130] In equation (3), h b denoted as riverbed elevation (m); x represents the distance along the river channel (m).

[0131] Friction slope S f This represents the energy loss caused by friction between the riverbed and the water body, which can generally be calculated using the Manning formula, expressed as:

[0132]

[0133] In equation (4), n is the Manning roughness (s / m) 1 / 3 ); It is the hydraulic radius (m), where P is the wetted perimeter (m).

[0134] The formula for calculating hydrostatic pressure I1 is:

[0135]

[0136] In equation (5), I2 is the wall pressure of the natural river channel, expressed as:

[0137]

[0138] In equation (6), h(x,t) is the water depth (m), which varies with time and space. η is the water level (m), and B(x,η) is the cross-sectional width (m) of the wetted perimeter of the open channel when the water level of the natural river is η.

[0139] The partial differential equation in its conserved form can be expressed as a governing equation, as follows:

[0140]

[0141] In equation (7), U is a vector of conserved variables. F is the flux vector. S is the source term vector.

[0142] The Jacobian matrix of the governing equations is expressed as:

[0143]

[0144] In equation (8), c is the gravitational wave velocity (m / s), u is the average flow velocity (m / s), and T is the width of the water surface (m).

[0145] Its eigenvalues ​​and eigenvectors are represented as follows:

[0146]

[0147] The eigenvalues ​​of the Jacobian matrices mentioned above are all real numbers and are two different values. Therefore, the governing equations with unequal real eigenvalues ​​constitute a hyperbolic system. In the problem domain of the calculation, discontinuous numerical solutions may appear. Therefore, the RKDG method uses the discontinuous finite element method to solve related problems.

[0148] Step 203: Based on the discontinuous finite element method and the Runge-Kutta time-discrete shallow water equations.

[0149] like Figure 3 As shown, the discontinuous finite element method divides the one-dimensional problem domain under study into Ne elements with Ne+1 nodes. It is assumed that the variables inside the elements are continuous, while the elements and their boundaries have discontinuous values.

[0150] The integral formula for the equation applied to each element in the discontinuous Galerkin method is as follows:

[0151]

[0152] In equation (10), and These are the coordinates of the starting and ending points of the finite element element, N. i and N j U is a shape function. j Let F represent the conserved variables (flow rate Q and cross-sectional area A), F represent the flux vector, and S represent the source term vector. This represents numerical flux, where N is a shape function.

[0153] The mass and momentum equations of the shallow water equations can both be written in discontinuous Galerkin form, as follows:

[0154]

[0155] Where z represents the water level.

[0156] The final form of the discontinuous finite element discrete one-dimensional shallow water equation is:

[0157]

[0158] In equation (13), Δx represents the length of the element, and ξ represents the local coordinates, which means mapping the physical element to ξ∈[-1,1], making it convenient to perform numerical integration using the Gaussian integral formula. i (ξ) represents a shape function. S i (ξ) represents the source term.

[0159] In the formula, p is the mass matrix, expressed as:

[0160]

[0161] In equation (14), Ω e Representing unit e, dΩ represents the integral over the e-th unit. dΩ indicates that the variable of integration is Ω, which is a quantity representing length. dΩ is an infinitesimal.

[0162] In terms of time, the RKDG method uses the third-order Runge-Kutta method (TVD Runge-Kutta) with reduced total variation for discretization, as shown in the following formula:

[0163]

[0164] In equation (15), Δt=min(Δt i ), where Δt represents the time step, Δt i This represents the time step of the i-th unit. CFL is an abbreviation for the Courant-Friedrichs-Lewy condition. It is an important stability condition used to ensure the computational stability of the numerical solution. n represents the n-th step in the time-progression calculation; its superscript indicates that the value is taken at the n-th time step.

[0165] Step 204: Solve for numerical flux by combining the inherently non-oscillating ENO slope limiter and the HLL Riemann solver.

[0166] To address the discontinuous numerical solutions and shock waves encountered in numerical calculations, this model employs an ENO-type slope limiter. The fundamental idea is to reconstruct variables with finite slopes using adjacent elements. The model assumes that the mean of variable U is... It is actually a first-degree polynomial, so it needs to be restructured. The other form of the polynomial is:

[0167] A = ax + b, Q = cx + d

[0168] Solving the equations of this unit and the adjacent units before and after it simultaneously yields the following results:

[0169]

[0170]

[0171] Solve the system of equations by simultaneously solving for a and b, c and d, and determine the slope and intercept of the linear polynomial with the smaller slope. e(i) represents the i-th unit. This represents the average cross-sectional area of ​​the water passage in the i-th unit. This represents the average flow rate in the i-th unit.

[0172] In the discontinuous Galerkin method of this model, the numerical flux at the element-to-element boundary is... The HLL flux function is used instead, and is expressed as:

[0173] Among them, F L F represents the flux of the cell on the left side of the interface. R This is represented as the flux of the right-hand cell.

[0174] In the formula, S L The wave velocity on the left side at the boundary between elements is expressed as:

[0175]

[0176] S R The wave velocity on the right side at the element-to-element boundary is expressed as:

[0177]

[0178] In the above formula, u* and c* are represented as follows:

[0179]

[0180] In the formula, L and R represent the limit values ​​of the parameter on both sides of the element boundary. - u+ represents the flow velocity to the left of the element boundary, and u+ represents the flow velocity to the right of the element boundary; h + h represents the water depth on the right side of the cell. - The value represents the water depth on the left side of the cell; u* represents the intermediate computational cost with the same flow velocity at the boundary; and c* represents the intermediate computational cost related to the wave velocity at the boundary.

[0181] The area calculation of each cross-section in this model uses the calculus dispersion theorem, which is highly accurate and easy to implement in programming. The specific integration formula is as follows:

[0182]

[0183] In the formula, x and y are the coordinates of the cross-section nodes, the subscripts 1 and 2 represent the left and right node numbers at any position, and n is the number of elements.

[0184] Furthermore, in this embodiment, the construction of the one-dimensional river hydrodynamic numerical model, through meticulous model parameter calibration, ensures that the model can highly reflect actual hydrological conditions, thereby achieving high-precision simulation of flow calculations. Based on this, the application boundaries of the model are expanded to explore forward and reverse applications of the hydrodynamic model, achieving bidirectional flow calculations from upstream to downstream and from downstream to upstream. This verifies the model's accuracy and generalization capabilities in calculations across different directions, providing a more comprehensive perspective on hydrodynamics.

[0185] Step 205: Perform forward and reverse verification on the above one-dimensional river hydrodynamic numerical model.

[0186] The positive verification mentioned above compares the simulated water depth of the one-dimensional river hydrodynamic numerical model with the theoretical water depth, and calculates the mean absolute error (MAE) and mean relative error (MRE) to evaluate the accuracy of the model.

[0187] For example, a forward simulation using a one-dimensional river hydrodynamic numerical model is performed, with the upper boundary being a constant flow rate of 157.9 m³. 3 / s, free outflow at the lower boundary, yielding water depths at each cross-section as follows: Figure 4 As shown, the simulated water depth is compared with the theoretical water depth using a plot. Mean absolute error (MAE) and mean relative error (MRE) are selected as analytical indicators to evaluate the accuracy of the hydrodynamic simulation, expressed as follows:

[0188]

[0189] The initial calculation yielded a mean absolute error (MAE) of 0.698. The mean relative error (MRE) was 22.7%, which is relatively large, necessitating model parameter calibration.

[0190] After model parameter calibration, simulation was performed again using a one-dimensional river hydrodynamic numerical model, yielding results at 157.9m. 3 The water depth and flow rate curves of each cross-section flowing uniformly down / s, such as Figure 5 , Figure 6 As shown, the simulated water depth was compared with the theoretical water depth, and the mean absolute error (MAE) was calculated to be 0.104. The mean relative error (MRE) was 3.4%. The results show that the small errors indicate that the model simulation results can reflect the actual water depth and are reliable.

[0191] The aforementioned reverse verification involves calculating the water flow from downstream to upstream, comparing the simulated water depth of the one-dimensional river hydrodynamic numerical model with the theoretical water depth, and calculating the mean absolute error (MAE) and mean relative error (MRE) to evaluate the model's accuracy.

[0192] For example, a one-dimensional river hydrodynamic numerical model is used for reverse simulation, that is, the water flow is calculated from downstream to upstream. The upper boundary is a constant water level of 8m (water depth of 3m), and the lower boundary is a constant flow rate of 157.9m. 3 / s, to obtain the water depth at each cross-section, such as Figure 7 As shown in the figure, the plotted water depths from the reverse simulation were compared with the theoretical water depths, and the mean absolute error (MAE) was calculated to be 0.251. The mean relative error (MRE) was 8.1%, which is relatively small, indicating that the reverse simulation can also reflect the actual water depth. However, the accuracy of the reverse simulation is not as good as that of the forward simulation. Therefore, the subsequent calculations used the forward simulation of the one-dimensional river hydrodynamic numerical model.

[0193] Step 106: Assume an initial tributary flow rate, and gradually adjust the tributary flow rate to make the simulated outlet flow rate process approximate the actual outlet flow rate process.

[0194] Specifically, assuming zero tributary flow, and considering only the known inflow process—the upstream "reservoir discharge"—a "tributary-free outflow" process is simulated using a one-dimensional river hydrodynamic model. The simulated "tributary-free outflow" is compared with the actual "station flow measurement," and the difference between the two is calculated. This calculated difference is considered the flow that the tributary should replenish. Based on the river's hydraulic characteristics and water propagation velocity, this flow difference is replenished to the tributary at an earlier time to correct the tributary's contribution to the outflow.

[0195] Through the above adjustments, the tributary flow rate is updated, and the simulation and comparison are performed again. This process iterates continuously, gradually reducing the error between the simulated and actual outlet flow rates until the preset accuracy is achieved. The core of the iterative method lies in using actual flow data to correct the tributary flow rate assumptions, thereby improving the model's accuracy. By iteratively correcting the tributary flow rate, the error between the simulated outlet flow rate and the actual observed flow rate is kept within the preset accuracy range. In each iteration, the tributary flow rate correction strategy is adjusted according to the tributary's inflow location to adapt to flow changes under different scenarios and dynamically optimize the simulation results.

[0196] Furthermore, to stabilize the model, the simulation was run for 200 hours based on the initial flow rate, building upon the initial 60 hours of data. Therefore, the upstream inlet flow rate for the first 200 hours was 157.9 m³ / s. 3 / s, "Reservoir discharge" data starts from hour 201. The relationship between "Reservoir discharge" at the upstream inlet, "Station flow measurement" at the downstream outlet, and the simulated "outflow without tributaries" is shown in [reference needed]. Figure 8 As shown.

[0197] After reservoir discharge and tributary flow into the river channel, the water requires a certain propagation time to reach the monitoring section before it can be observed at the final section, resulting in a time lag response problem. Figure 9 As shown. The outflow value Q at section B. B (t B The inflow value Q at section A A (t A The following time lag relationship exists between them, expressed as:

[0198] t B =t A +T A-B (29)

[0199] In the formula, T A-B T is the time (i.e., the delay) for the flow rate at section A to propagate to section B. A-B The magnitude can be determined by applying a certain disturbance to the inflow value at section A and then observing the time when the maximum fluctuation occurs at section B.

[0200] Assuming the reservoir discharge is uniform flow, the lag response time of the monitoring section is determined by artificially applying different degrees of disturbance to the inlet section (increasing the flow rate by different amounts at time 0). The disturbance applied to the inlet section's flow rate during the monitoring section's response process is as follows: Figure 9 As shown. From Figure 10 As can be seen, although the degree of disturbance varies, the peak time displayed by the monitoring section is around 22 hours, at which point the flow rate observed at the monitoring section contains the most inflow information. Therefore, the lag time corresponding to the reservoir discharge is taken as 22 hours.

[0201] from Figure 8 It can also be seen that the peak time of the inlet flow is at 226h, while the peak time of the "outflow without tributaries" process generated by the flow propagating to the monitoring section is at 249h. This verifies that it takes about 22 hours for the reservoir discharge to propagate to the monitoring section. The lag time of the reservoir discharge can be taken as 22h, and this generalization method is reliable.

[0202] In this embodiment, a one-dimensional river hydrodynamic model is constructed to calculate the tributary confluence flow at different locations within the river channel based on the flow process observed at the river channel flood control section and the reservoir discharge process. The accuracy of the model simulation results is evaluated based on the difference between the simulated outflow and the station-measured flow. Specifically, the river channel is divided into multiple calculation intervals according to the different tributary confluence locations. Within each calculation interval, the flow changes before and after the tributary confluence are simulated using the one-dimensional hydrodynamic model. For each calculation interval, the spatiotemporal distribution characteristics of the interval flow are calculated based on the simulation results. The relationship between the interval flow and the reservoir discharge flow is evaluated by combining the flow changes in different calculation intervals. A disturbance is applied at the tributary confluence point, and the tributary lag time is determined by the response time of the monitoring section. The input time of the tributary flow is adjusted according to the lag time to simulate the impact of the tributary on the river channel flow.

[0203] The following section describes the interval flow calculation method provided in this application under three scenarios: tributaries flowing in from 5km upstream, 50km midstream, and 95km downstream.

[0204] When a tributary joins the river 5 km upstream, and the river channel is known to be 100 km long, assuming the reservoir discharge is uniform flow, the lag response time of the monitoring section is determined by artificially applying different degrees of disturbance to the tributary cross-section (increasing the flow rate by different amounts at time 0). The response process of the disturbance applied to the tributary cross-section at the monitoring section is as follows: Figure 11 As shown. From Figure 11 As can be seen, although the degree of disturbance varies, the peak time displayed by the monitoring section is around 22 hours. At this time, the flow value observed by the monitoring section contains the most information about the tributary confluence. Therefore, the lag time corresponding to this tributary is taken as 22 hours.

[0205] First, assume the tributary flow is zero. Since the reservoir discharge delay is 22 hours, start calculating the difference between "no tributary outflow" and "station flow measurement" from hour 223. Because the tributary delay is 22 hours, advance the calculated difference by 22 hours, starting tributary replenishment from hour 201. Figure 12 As shown. Then, the new tributary flow is added to the model, the simulated new outflow process is subtracted from the "site flow measurement", and the difference is replenished to the tributary starting from 201h, and so on.

[0206] Analysis revealed that the initial iterations of the hydrodynamic model for tributary flow showed significant errors between the simulated and actual outlet flow rates, with both peak values ​​and peak occurrence times deviating from the measured flow patterns. By the third iteration, the mean relative error (MRE) between the simulated and actual outlet flow rates converged to within 1%, indicating a marked improvement in flow deviation. To further reduce the error and make the simulated outflow more closely resemble actual station measurements, six iterations were performed, reducing the maximum difference between the simulated and actual outlet flow rates to 3.8 m. 3 / s, the mean absolute error (MAE) and mean relative error (MRE) decreased to 1.51% and 0.39%, respectively. Figure 13 As shown. At this point, the tributary joins the river channel from a point x = 5 km upstream, and the tributary flow process is as follows. Figure 14 As shown.

[0207] When a tributary joins the river 50 km downstream, and the river channel is known to be 100 km long, assuming uniform flow from the reservoir, the lag response time of the monitoring section is determined by artificially applying different degrees of disturbance to the tributary cross-section (increasing the flow rate by different amounts at time 0). The response process of the disturbance applied to the tributary cross-section at the monitoring section is as follows: Figure 15 As shown. From Figure 15 As can be seen, although the degree of disturbance varies, the peak time displayed by the monitoring section is around 8 hours. At this time, the flow value observed by the monitoring section contains the most information about the tributary confluence. Therefore, the lag time corresponding to this tributary is taken as 8 hours.

[0208] First, assume the tributary flow is zero. Since the reservoir discharge delay is 22 hours, start calculating the difference between "no tributary outflow" and "station flow measurement" from hour 223. Because the tributary delay is 8 hours, advance the calculated difference by 8 hours, starting tributary replenishment from hour 215. Figure 16 As shown. Then, the new tributary flow is added to the model, and the simulated new outflow process is subtracted from the "site flow measurement". Similarly, the difference is added to the tributary starting from 215h, and so on.

[0209] Analysis revealed that the initial iterations of the hydrodynamic model for tributary flow showed significant errors between the simulated and actual outlet flow rates, with both peak values ​​and peak times deviating from the measured flow patterns. By the third iteration, the mean relative error (MRE) between the simulated and actual outlet flow rates converged to within 1%, indicating a marked improvement in flow deviation. To further reduce the error and make the simulated outflow more closely resemble actual station measurements, five iterations were performed, reducing the maximum difference between the simulated and actual outlet flow rates to 3.3 m. 3 / s, the mean absolute error (MAE) and mean relative error (MRE) decreased to 1.29% and 0.34%, respectively. Figure 17As shown in Figure 18. At this point, the tributary joins the river channel at a point x = 50 km upstream, and the tributary flow process is shown in Figure 18.

[0210] When a tributary joins the river 95 km downstream, and the river channel is known to be 100 km long, assuming the reservoir discharge is uniform, the lag response time of the monitoring section is determined by artificially applying different degrees of disturbance to the tributary cross-section (increasing the flow rate by different amounts at time 0). The response process of the disturbance applied to the tributary cross-section at the monitoring section is as follows: Figure 19 As shown. From Figure 19 As can be seen, although the degree of disturbance varies, the peak time displayed by the monitoring section is around 0h. At this time, the flow value observed by the monitoring section contains the most tributary confluence information. Therefore, the lag time corresponding to this tributary is taken as 0h.

[0211] First, assume the tributary flow is zero. Since the reservoir discharge delay is 22 hours, start calculating the difference between "no tributary outflow" and "station flow measurement" from hour 223. Because the tributary delay is 0 hours, directly replenish the tributary with the calculated difference starting from hour 223. Figure 20 As shown. Then, the new tributary flow is added to the model, the simulated new outflow process is subtracted from the "site flow measurement", and the difference is replenished to the tributary starting from 223h, and so on.

[0212] Analysis revealed that the simulated tributary flow, obtained only through initial iterations of the hydrodynamic model, exhibited significant errors compared to the actual outlet flow, with both peak value and peak occurrence time deviating from the measured flow rate. However, by the third iteration, the mean relative error (MRE) between the simulated and actual outlet flow converged to within 1%, indicating a marked improvement in flow deviation. To further reduce the error and make the simulated outflow more closely resemble the actual flow measurements at the site, four more iterations were performed, reducing the maximum difference between the simulated and actual outlet flow to 3.0 m. 3 / s, the mean absolute error (MAE) and mean relative error (MRE) decreased to 0.83% and 0.22%, respectively. Figure 21 As shown. At this point, the tributary joins the river channel from a downstream location x = 95 km away, and the tributary flow process is as follows. Figure 22 As shown.

[0213] This application's embodiments calculate the interval flow characteristics, including peak flow, peak occurrence time, and average flow, based on the different tributary confluence locations at upstream, midstream, and downstream levels. By comparing the interval flow characteristics under different scenarios, the impact of tributary confluence location on interval flow is analyzed. Figure 23The diagram illustrates the flow patterns of the tributary confluence under three scenarios: 5km upstream, 50km midstream, and 95km downstream. It can be seen that as the tributary's confluence location changes, the closer to the downstream outlet, the smaller the tributary's flood peak and the later its occurrence. Figure 24 As shown.

Claims

1. A method for calculating river flow, characterized in that, include: Acquire topographic, hydrological, and roughness data to construct a one-dimensional river channel hydrodynamic numerical model; The one-dimensional river hydrodynamic numerical model is based on the discontinuous finite element method and the Runge-Kutta time-discrete shallow water equations, and combines the inherently non-oscillating ENO slope limiter and the HLL Riemann solver to solve the numerical flux. For different tributary confluence scenarios, the interval flow is calculated using an iterative method. The iterative method includes an initial assumption that the tributary flow is zero, and the simulation of tributary-free outflow using the one-dimensional river hydrodynamic numerical model is used to determine the difference between the simulated tributary-free outflow and the actual flow measurement at the station. The tributary flow rate is adjusted based on the differences, river hydraulic characteristics, and water propagation speed, and the process is repeated iteratively until the error between the simulated outlet flow rate and the actual outlet flow rate reaches the preset accuracy.

2. The method for calculating river flow according to claim 1, characterized in that, The governing equations of the model are the Saint-Venant equations for one-dimensional shallow water flow, including the continuity equation and the motion equation, which are used to describe the continuity of river flow and the conservation of momentum. The continuity equation is expressed as: In the formula, A is the cross-sectional area of ​​the water passage; Q is the cross-sectional flow rate; The equation of motion is expressed as: In the formula, S x S represents the riverbed slope. f I1 is the friction slope; I2 is the hydrostatic pressure; g is the wall pressure of the natural river channel; g is the acceleration due to gravity. Among them, the riverbed slope S x The calculation formula is: In the formula, h b denoted as the riverbed elevation, and x as the distance along the river channel. Friction slope S f This represents the energy loss caused by friction between the riverbed and the water body, calculated using Manning's formula, and is expressed as: In the formula, n is the Manning roughness. It is the hydraulic radius, where P is the wetted perimeter; The formula for calculating hydrostatic pressure I1 is: The formula for calculating the wall pressure I2 of a natural river channel is: Where h is the water depth, which varies with time and space; η is the water level; B(x,η) is the cross-sectional width of the wetted perimeter of the open channel when the water level of the natural river is η. The partial differential equation in its conserved form can be expressed as a governing equation as follows: In equation (7), U is a vector of conserved variables. F is the flux vector. S is the source term vector. The Jacobian matrix of the governing equations is expressed as: In the formula, c is the gravitational wave velocity; u is the average flow velocity; and T is the width of the water surface. Its eigenvalues ​​and eigenvectors are represented as follows:

3. The method for calculating river flow according to claim 1, characterized in that, After constructing the one-dimensional river hydrodynamic numerical model, the method further includes performing forward and reverse verification on the one-dimensional river hydrodynamic numerical model. The positive verification involves comparing the simulated water depth of the one-dimensional river hydrodynamic numerical model with the theoretical water depth, and calculating the mean absolute error (MAE) and mean relative error (MRE) to evaluate the accuracy of the model. The reverse verification is performed by calculating the water flow from downstream to upstream, comparing the simulated water depth of the one-dimensional river hydrodynamic numerical model with the theoretical water depth, and calculating the mean absolute error (MAE) and mean relative error (MRE) to evaluate the accuracy of the model.

4. The method for calculating river flow according to claim 1, characterized in that, The one-dimensional river hydrodynamic numerical model is based on the discontinuous finite element method and the Runge-Kutta time-discrete shallow water equations, including: The one-dimensional problem domain is divided into Ne elements and Ne+1 nodes using the discontinuous finite method. Variables within an element are continuous, while those at element boundaries are discontinuous. The integral formula for the equation applied to each element in the discontinuous Galerkin method is as follows: In the formula, and These are the coordinates of the starting and ending points of the finite element element, N. i and N j U is a shape function. j Let F represent the conserved variables of flow rate Q and cross-sectional area A, F represent the flux vector, and S represent the source term vector. Represents numerical flux; N is a shape function; The mass and momentum equations of the shallow water equations are both converted to discontinuous Galerkin form, expressed as: Where z represents the water level; The discontinuous finite element method for discretizing the one-dimensional shallow water equations is expressed as: In the formula, Δx represents the length of the element, ξ represents the local coordinates, that is, mapping the physical element to ξ∈[-1,1], which facilitates numerical integration using the Gaussian integral formula; N i (ξ) represents a shape function; S i (ξ) represents the source term; In the formula, p is the mass matrix, expressed as: In the formula, Ω e Representing unit e, This represents the integral over the e-th unit; Discretization in time is performed using the Runge-Kutta method, which reduces the total variation by third order, and is expressed as follows: In the formula, Δt represents the time step, Δt i The value represents the time step of the i-th unit. CFL stands for Courant-Friedrichs-Lewy condition, which is used to ensure the computational stability of the numerical solution. n represents the n-th step in the time-progression calculation, and the superscript indicates that the value is taken at the n-th time step.

5. The method for calculating river flow according to claim 1, characterized in that, Solving numerical fluxes by combining the inherently non-oscillating ENO slope limiter and the HLL Riemann solver includes: By reconstructing variables with finite slopes using adjacent elements, the equations of this element and its preceding and following adjacent elements are solved simultaneously to obtain: By solving the system of equations simultaneously for a and b, c and d, and determining the slope and intercept of the first-order polynomial with the smaller slope, e(i) represents the i-th unit. This represents the average cross-sectional area of ​​the water passage in the i-th unit. This represents the average flow rate in the i-th unit; In the Galerkin method for discontinuous hydrodynamics in one-dimensional channels, the numerical flux at the boundary between elements is... The HLL flux function is used instead, and is expressed as: Among them, F L F represents the flux of the cell on the left side of the interface. R This is represented as the flux of the right-hand cell; In the formula, S L The wave velocity on the left side at the boundary between elements is expressed as: S R The wave velocity on the right side at the element-to-element boundary is expressed as: Where u* and c* are respectively represented as: In the formula, L and R represent the limit values ​​on both sides of the element boundary, respectively; u - u+ represents the flow velocity to the left of the element boundary, and u+ represents the flow velocity to the right of the element boundary; h + h represents the water depth on the right side of the cell. - The value represents the water depth on the left side of the cell; u* represents the intermediate computational cost with the same flow velocity at the boundary; and c* represents the intermediate computational cost related to the wave velocity at the boundary.

6. The method for calculating river flow according to claim 1, characterized in that, The calculation of the interval flow includes: The river channel is divided into multiple calculation intervals based on the location where the tributaries flow into it; Within each calculation interval, the flow changes before and after the tributary's inflow are simulated using the one-dimensional hydrodynamic model. For each calculation interval, the spatiotemporal distribution characteristics of the interval flow are calculated based on the simulation results; By combining the flow changes in different calculation intervals, the relationship between interval flow and reservoir discharge flow is evaluated; A disturbance is applied at the confluence of the tributaries, and the tributary lag time is determined by monitoring the response time of the cross section. The input time of the tributary flow is then adjusted based on the lag time to simulate the impact of the tributaries on the river flow.

7. The method for calculating river flow according to claim 6, characterized in that, Depending on the location where the tributary flows into the area, the calculation of the interval flow also includes: Calculate the inter-regional flow characteristics at the upstream, midstream and downstream confluence points, including peak flow, peak occurrence time and average flow. By comparing the characteristics of inter-regional flow under different scenarios, the influence of the location of tributary confluence on inter-regional flow is analyzed.

8. The method for calculating river flow according to claim 1, characterized in that, Repeat the iteration until the error between the simulated outlet flow and the actual outlet flow reaches the preset accuracy, including: The tributary flow is gradually corrected by iterative method so that the error between the simulated outlet flow and the actual observed flow is within the preset accuracy range. In each iteration, a correction strategy is adopted to adjust the tributary flow based on the tributary inflow location, in order to adapt to the water flow changes under different scenarios and dynamically optimize the simulation results.

Images