The present application claims priority benefit to Chinese application No. 202111429722.5 filed on Nov. 29, 2021, the content of which is hereby incorporated by reference in its entirety.
The invention relates to the field of hydraulic engineering, in particular to the field of flood control risk analysis, and more particularly, to a numerical simulation method for the overflowing of river channel structures based on finite volume method.
A river channel structure refers to a structure that exists in a river channel, such as retaining weirs, river sluices, bridges and culverts, which have a significant impact on the water movement of the river. In engineering practice, a one-dimensional mathematical model is usually used to simulate the flow propagation occurred in the rivers. However, the water flow regime at the structure is extremely complex, and the control equation set used in the one-dimensional model cannot describe the complex water flow regime accurately. The structure is also known as a singularity within the computational domain and needs to be processed dedicatedly. In engineering practice, structures in river channels are very common. Therefore, whether the one-dimensional model can reasonably simulate the overflowing process of structures in the river channel determines whether the model can be popularized and applied on a large scale.
In the existing one-dimensional numerical simulation methods, the implicit finite difference method is regarded as the most popular numerical method for river flood simulation as it is developed earlier, originated in the 1960s. Common finite difference formats include Preissmann scheme, Abbott-Lonescu scheme, etc. Dealing with the structures in such model is completed by adding virtual sections and solving large sparse matrices. However, such numerical schemes based on implicit finite difference have no ability to capture shock waves and process excessive flow regimes. In some river channels with steep terrains and large slope changes, the numerical solution often diverges and even the simulation fails.
In recent years, the Godunov finite volume method based on solving the Riemann approximate solution has been widely applied in the field of flood numerical simulation. It can not only simulate the flow of river with a gentle slope, but also adaptively simulate the change of flow regime, and can also be well adapted to the river channel in a hilly area with a steep slope. Therefore, such method has been gradually developed into a new calculation scheme for simulating river flow. How to deal with the overflowing problem of river channel structures under the numerical framework of the new finite volume method is also a hot and difficult problem in this field at present.
The objective of the invention is to provide a numerical simulation method for the overflowing of river channel structures based on finite volume method. According to this method, the overflowing characteristics of all structures can be accurately reflected under the numerical framework of the one-dimensional explicit finite volume method, while the precision and stability of overall calculation of the one-dimensional river channel flow model are guaranteed.
The invention is implemented through the following technical solutions.
A numerical simulation method for the overflowing of river channel structures based on finite volume method, which uses a one-dimensional finite volume element to discrete the river channel, and takes the position of the river channel structures as the special interface of the river calculation element. The amount of water passing through the interface in each time step is calculated according to the structure overflowing formula (e.g., weir flow formula, gate outflow formula, etc.). The amount of water is then respectively deducted from and increased to an upstream element and a downstream element adjacent to the interface by adopting a source item processing method. The primitive variable values at the special interface are numerically reconstructed in order to ensure the continuity of calculation. The method includes the following specific steps:
1) Obtaining basic data of river channels and structures: obtaining the plane geometric shape and control section shape data of the river channel, wherein the distance between adjacent sections is not greater than 1 km; obtaining the spatial position information and geometric size information of the river channel structure, and arranging two control sections on the upstream and downstream of the structure, wherein a distance between the sections is not more than 100 m;
2) Performing spatial discretization of the river channel: one-dimensional finite volume element is used to discretize the river channel, wherein the position of each section is the center position of the element, and the position of the midpoint of two sections is the element interface position. If there is no structure at the position, the interface is referred to as a conventional interface; if there is a structure at the position, the interface is referred to as a special interface. The primitive variable values of the river are stored in the center of the elements;
3) Initializing the calculation conditions: assigning initial primitive variable values to each element of the river channel, i.e., initial water level and initial flow value, setting the initial operating state of the structures, and obtaining the initial values of the upstream and downstream boundaries of the river channel;
4) Obtaining the outer boundary conditions at time t, and obtaining the calculation time step dt according to the CFL (Courant-Friedrichs-Lewy) condition;
5) Solving the numerical flux at each conventional interface, the overflow of the structure, and the numerical flux at the special interface: describing the river channel flow movement by using the Saint-Venant equations with source term; calculating the numerical flux at the conventional interface of each element at time t by using the finite volume method based on the HLL (Harten-Lax-van Leer) approximate Rieman solution; calculating the value of flow passing through the structure at time t by applying the primitive variable values of the upstream and downstream elements of the structure according to overflowing characteristics of structures (the overflowing characteristics are, for example, as follows: some structures are only controlled by the upstream water level, some are controlled by the water levels in upstream and downstream at the same time, some of the water flow passing through the structure are open channel flow, some are pressure flow, etc.), and processing the flow value as a source term of a continuity equation in the upstream and downstream elements; reconstructing the primitive variable values at the structure interface by using a non-reflection boundary condition in order to calculate the continuity of calculation, and then calculating the numerical flux passing through the special interface;
6) Obtaining the primitive variable values of each element at time t+dt: updating the primitive variable values at the center of each element at time t+dt through a numerical flux value at each element interface at time t and a source term value of overflowing of structures, and updating boundary conditions at both ends of the river channel synchronously to time t+dt;
7) Setting t=t+dt, and repeating the steps 4) to 6) until the end of the calculation. Further, in the step 4), the selection of dt is limited by the CFL condition, as shown in Formula (1):
in which, Ncfl is the CFL number, u is the average flow velocity of the section, c is the wave velocity, Δx is the space step of the finite volume element, and dt is the time step.
Further, in the step 5), one-dimensional unsteady flows in a river channel considering the overflowing of the structure are described by the Saint-Venant equations with source term, as shown in Formula (2):
in which, x is the space variable, t is the time variable, and D, U, F, and S are vector representations of variables in the equation set, to be specific:
in which, B is the width of the water surface, Z is the water level, Q is the discharge Z and Q are called primitive variables, A is cross-sectional area, f1 and f2 represent two components of the vector F(U), respectively, g is the acceleration of gravity, J is the on-way resistance loss, with expression of J=(n2Q|Q|)/A2R4/3), R is the hydraulic radius, n is the Maiming roughness coefficient, and q1 is the source term value per unit length of the river channel.
Further, in the step 5), the overflow through the structure at time t is calculated through the primitive variable values of the upstream and downstream elements of the structure according to the overflowing characteristics of the structure. Taking the retaining weir in the river channel as an example, the overflowing calculation formula is shown in Formula (3):
in which, hup=Zup−Zweir; hdown=Zdown−Zwein; qw is the flow discharge passing through the structure; Zup and Zdown are the water levels of calculation elements on the upstream and downstream of the structure, respectively; Zweir is the weir crest elevation; and l is the weir width.
Further, in the step 5), the primitive variable values at the special interface are calculated by using the non-reflection boundary condition in order to ensure the continuity of calculation of all interface fluxes, to be specific: when the numerical fluxes of the upstream and downstream elements of the structure passing through the interface are calculated respectively, the primitive variable values at the interface are reconstructed as:
Z*=Z
up
,Q*=0 (4)
Z*=Z
down
,Q*=0 (5)
in which, are the water level and the discharge after reconstruction at the element interface, respectively.
The invention has the following advantages and beneficial effects.
The method proposed in the invention can reflect the overflowing characteristics of each structure. Meanwhile, the precision and stability of the one-dimensional river channel flow model can also be guaranteed, thereby providing a new solution for dealing with overflowing process of structures under the numerical framework of a one-dimensional finite volume method.
The invention will be further described in conjunction with the attached figures and embodiments.
in which, i−1, i and i+1 are the numbers of finite volume elements; i−1/2 and i+1/2 are left and right interfaces of the ith element Δxi is the length of the ith finite volume element; and F*i−1/2 and F*i+1/2 are numerical fluxes passing through the left and right interfaces of the element;
The invention will be further described in conjunction with
The invention provides a numerical simulation method for the overflowing of river channel structures based on finite volume method, which uses a one-dimensional finite volume element to discrete the river channel, and takes the position of the river channel structures as the special interface of the river calculation element. The amount of water passing through the special interface in each time step is calculated according to the structure overflowing formula (e.g., weir flow formula, sluice gate outflow formula, etc.). The amount of water is then respectively deducted from and increased to an upstream element and a downstream element adjacent to the interface by adopting a source item processing method. The primitive variable values at the special interface are numerically reconstructed in order to ensure the continuity of calculation. The method can accurately reflect the overflowing characteristics of all structures under the numerical framework of the one-dimensional explicit finite volume method, while the precision and stability of overall calculation of the one-dimensional river channel flow model are guaranteed. The method includes the following specific steps:
1) Obtaining the plane geometric shape and control section shape data of the river channel, wherein the distance between adjacent sections is not greater than 1 km; obtaining the spatial position information and geometric size information of the river channel structure, and arranging two control sections on the upstream and downstream of the structure, wherein a distance between the sections is not more than 100 m;
2) One-dimensional finite volume element is used to discretize the river channel, wherein the position of each section is the center position of the element, and the position of the midpoint of two sections is the element interface position. If there is no structure at the position, the interface is referred to as a conventional interface; if there is a structure at the position, the interface is referred to as a special interface. The primitive variable values of the river are stored in the center of the element. The details are plotted in
3) Initializing the calculation conditions, assigning initial primitive variable values to each element of the river channel, i.e., initial water level and initial discharge, setting the initial operating state of the structures, and obtaining the initial values of the upstream and downstream boundaries of the river channel;
4) Obtaining the calculation time step dt according to the CFL (Courant-Friedrichs-Lewy) condition. The specific conditions of CFL are shown in Formula (1):
in which, u is the average flow velocity of the section, c is a wave velocity, Δx is the space step of the finite volume element, and dt is the time step. In order to ensure the stability of the overall numerical calculation, Ncfl is recommended to be not exceeding 1.0.
5) Describing one-dimensional flows movement in a river channel considering the overflowing of the structure by using the Saint-Venant equations with source term, as shown in Formula (2):
in which, x is the space variable, t is the time variable, and D, U, F, and S are vector representations of variables in the equation set, to be specific:
in which, B is the width of the water surface, Z is the water level, Q is the discharge, Z and Q are called primitive variables, A is the cross-sectional area, f1 and f2 represent two components of the vector F(U), respectively, g is the acceleration of gravity, J is the on-way resistance loss, with expression of J=(n2Q|Q|)/(A2R4/3), R is the hydraulic radius, n is the Manning roughness coefficient, and q1 is the source term value per unit length of the river channel.
The numerical flux at the conventional interface of each element at time t is calculated by the finite volume method based on the HLL approximate Riemann solution. The detailed solution process of the HLL approximate Riemann solution can be found in the literature (Zhang Dawei, Cheng Xiaotao, Huang Jinchi, etc., “Widely adaptable numerical model for complicated open channel flows”, [J]. Journal of Hydraulic Engineering, 2010,41(4):531-536 (in Chinese), the content of which is herein incorporated by reference in its entirety; Zhang Dawei, Numerical Simulation of Dam Burst Flow Based on Godunov Format [M]. China Water&Power Press, Beijing, 2014,12, the content of which is herein incorporated by reference in its entirety).
The overflow through the structure at time t is calculated through the primitive variable values of the upstream and downstream elements of the structure according to the overflowing characteristics of the structure Taking the retaining weir in the river channel as an example, the overflowing calculation formula is shown in Formula (3):
in which, hup=Zup−Zweir; hdown=Zdown−Zweir; qw is the flow discharge passing through the structure; Zup and Zdown are the water levels of calculation elements on the upstream and downstream of the structures; Zweir is the weir crest elevation; and l is the weir width.
For different structures, the overflow at the current time can be obtained through the primitive variable values of their adjacent upstream and downstream elements. It should be noted that some structures have a definite overflowing formula can be used directly. While some structures are complicated in design, with no definite overflowing formula can be used. In this case, it is necessary to first acquire an overflowing curve formula of the structures through physical tests and other means, and then obtain the corresponding overflow based on the primitive variable values of the upstream and downstream of the structure. For a certain structure, the amount of water flowing from the upstream element of the structure at time t will inevitably enter the adjacent downstream element. The overflowing problem of the structure is considered by adding a source term in the governing equations, which can theoretically ensure the overall mass conservation characteristics of the mathematical model calculation. Meanwhile, this method can also ensure the stability of the calculation model.
The primitive variable values at the special interface are reconstructed by using the non-reflection boundary condition in order to ensure the continuity of calculation of all interface fluxes, to be specific: when the numerical fluxes of the upstream and downstream elements of the structure passing through the interface are calculated respectively, the primitive variable values at the interface are reconstructed as:
Z*=Z
up
,Q*=0 (4)
Z*=Z
down
,Q*=0 (5)
in which, Z*, Q* are the primitive variable values at the element interface, Zup is the water level of the upstream element of the structure, and Zdown is the water level of the downstream element of the structure.
According to the reconstructed primitive variable values at the special interface, the numerical flux through the special interface is further calculated.
6) Updating the primitive variable values at the center of each element at the time t+dt through a numerical flux value at each element interface at the time t and a source term value of overflowing of structures, and updating boundary conditions at both ends of the river channel synchronously to the time t+dt;
7) Setting t=t+dt, and repeating the steps 4) to 6) until the end of the calculation.
The above-mentioned embodiments are only a partial expression of the invention, and cannot cover the whole invention. On the basis of the above-mentioned embodiments and figures, those skilled in the field can obtain more embodiments without paying creative work. Therefore, these embodiments obtained without paying creative work should be included in the protection scope of the invention.
Number | Date | Country | Kind |
---|---|---|---|
202111429722.5 | Nov 2021 | CN | national |