Method for simulating braid-shaped river channel water flow plane flow field

Abstract

The invention discloses a method for simulating a braid-shaped river channel water flow plane flow field. The method comprises the following steps that a river channel two-dimensional calculation grid is manufactured; interpolating the topography of the grid points, and completing fine modeling of the plait-shaped river channel plane area riverbed topography and landform; carrying out numerical discretization on the two-dimensional hydrodynamic control equation on the two-dimensional computational grid, and constructing a basic two-dimensional hydrodynamic solver; a dry and wet dynamic boundary processing module in a memory and compensation mode is established, so that the dry and wet dynamic boundary processing module and a hydrodynamic solver work in a coordinated and synchronous mode, and a complete two-dimensional hydrodynamic model is jointly established; for a target river reach of the braided river channel, setting inlet and outlet flow and water level open boundary conditions, simulating a water flow plane flow field of the braided river channel, storing water level distribution and flow field results calculated by the model in a hard disk, and displaying the results on a computer screen for feedback to users. The method solves the problems that in the technical field of river simulation, a braided river flow plane flow field simulation method is too large in error, poor in effect or poor in method universality and the like.

Description

Technical Field

[0001] The present application relates to the field of hydraulic simulation, and specifically to a method for simulating the plane flow field of water flow in a braided river. Background Art

[0002] Common types of river channels in nature include straight, curved, bifurcated and braided. Braided rivers are widely distributed in the middle and lower reaches of river basins, deltas and other areas. They are the most complex type of river channels among all types of river channels, with wide and shallow multi-branched (interwoven like braids) riverbed characteristics and staggered and scattered water flow patterns. The difficulty of numerical simulation of the plane flow field of braided rivers is also the greatest among all types of rivers. At the same time, simulating and studying the plane flow field of braided river flow is also the basis for carrying out braided river regulation and river control, regional flood prevention and other work, and has important practical value in water conservancy engineering practice. The simulation of the plane flow field of the river channel generally adopts a planar two-dimensional hydrodynamic model (referred to as a two-dimensional model). The difficulty of using it to simulate the movement of water flow in braided rivers is analyzed as follows.

[0003] First, the width of braided channels is often more than 2~3km, and many of the grooves are only 50m wide or even smaller. The wide river means a large simulation area and more two-dimensional grid units are required. At the same time, accurate description of the distribution of shoals and narrow grooves requires the use of small-scale grids. Considering both the coverage of wide river areas and the accurate description of river regimes, the two-dimensional grid generation of braided channels (which belongs to the first part of two-dimensional modeling) has brought considerable difficulties. Moreover, in the conventional measured topographic maps of wide rivers (1:5000~1:10000), the spacing between elevation measurement points is usually more than 50~100m, which means that braided channels often have only a single or no topographic measurement point in the lateral direction of many of their grooves. Therefore, the spatial resolution of the conventional measured topography of braided channels is not sufficient to support the accurate interpolation of the two-dimensional grid topography of the channel (which belongs to the second part of two-dimensional modeling). These two factors make the two-dimensional fine modeling of braided channels very difficult. Second, the dry-wet state transition in braided channels is very frequent under the background of complex flow paths and changeable flow patterns, which puts forward stringent requirements on the simulation capability of the two-dimensional model for the wet-dry dynamic boundary. Third, the minimum size of the unit of the computational grid used to describe the complex braided channels is generally very small. Under this background, the alternating slow and rapid currents and scattered water flow in braided channels will increase the risk of instability in the calculation of the two-dimensional model, which puts forward high requirements on the numerical solution of the model.

[0004] Hydrodynamic models are divided into explicit and implicit types. Although the implicit hydrodynamic model has a complex structure, it has good numerical stability. With the rapid development of computational mathematics and computational fluid dynamics methods in this century, the semi-implicit hydrodynamic model based on the Euler-Lagrange method can be applied to simulate various slow, rapid and complex flow states and can use a large time step, which can better solve the problem of instability in the two-dimensional hydrodynamic calculation of braided rivers. However, at the same time, the problems of fine modeling of the two-dimensional model of braided rivers and the treatment of complex and frequent dry-wet dynamic boundaries have not been well solved, especially the lack of a dry-wet dynamic boundary treatment method with high stability and mass conservation that is suitable for the two-dimensional model of complex rivers.

[0005] The water area and its boundaries of the river channel will change frequently as the flood spreads, which is called the wet-dry dynamic boundary problem. In the calculation process of the two-dimensional hydrodynamic model, the change of the water area will cause some cells in the grid to be submerged when the water rises and dry out when the water falls. The calculation work used to realize the wet-dry conversion of the grid in order to characterize the change of the water area and its boundary is called wet-dry dynamic boundary simulation. The critical water depth method is the most commonly used method to simulate the wet-dry dynamic boundary. It defines a critical water depth h0 (for example, 0.001m) to define the wet-dry state of the grid element. When its water depth is less than h0, the grid element is dry, otherwise it is wet; at the same time, a wet-dry flag is defined for each grid element (one-to-one correspondence with the grid element) to record their current wet-dry state. The wet-dry state of the grid element is updated in real time according to the latest water flow information, so as to describe the dynamics of the water area and its boundary. In the traditional critical water depth method, when the water depth of a wet unit is less than h0, it is immediately "dried out" and exits the model calculation; when a dry unit meets the conditions for resuming flow, it is made wet and added to the model calculation at the next time step. The restoration flow conditions for the "dry→wet" conversion of the dry cell are as follows: ① The water level at the center of the dry cell is obtained by interpolating the information of the adjacent wet cells and the water depth is calculated, which should be greater than h0; ② There should be an inflow of water flux on the side of the dry cell, or there should be a large water level gradient between the dry cell and the adjacent wet cell.

[0006] The traditional critical water depth method cannot handle the potential negative water depth problem in the hydrodynamic calculation in time, which easily causes water volume non-conservation and distortion of calculation results under complex river conditions. This problem is particularly prominent in the implicit model with large time steps. The advantage of the implicit hydrodynamic model is that it allows the use of large time steps to carry out river flow calculations, but the grid wet and dry state is generally not updated in real time during the iterative solution of a single time step. These characteristics make it very easy for the implicit model based on the traditional critical water depth method to generate negative water depth units in the local area when encountering water areas with complex terrain and large water level gradients when using large time steps for calculations, which in turn causes water volume non-conservation and distortion of calculation results. Therefore, there is a contradiction between the large time step hydrodynamic calculation and the small time step required to simulate the dynamics of the rapidly changing water boundary. This problem has not been well solved so far and has caused large errors in the simulation of braided river flow. Summary of the invention

[0007] The purpose of the embodiments of the present application is to provide a method for simulating the plane flow field of water flow in a braided river, so as to solve the problems of excessive errors, poor results, and low versatility of the method in the current field of hydraulic simulation (river simulation) technology.

[0008] To achieve the above objectives, this application provides the following technical solutions:

[0009] The present invention provides a method for simulating a planar flow field of a braided river flow, comprising the following steps:

[0010] Step 1: Create a two-dimensional computational grid for the river channel, including classifying and marking the plane area of ​​the braided river channel, drawing the boundary line inside and outside the river channel, and generating computational grids by partition;

[0011] Step 2: interpolating the terrain of the grid points, including classifying and marking the scattered points and grid points of the braided river terrain, pairing the grid points with the scattered points of the terrain, and interpolating the terrain of the grid points, thereby completing the fine modeling of the riverbed topography and geomorphology in the planar area of ​​the braided river;

[0012] Step 3: numerically discretize the two-dimensional hydrodynamic control equations on a two-dimensional computational grid and construct a basic two-dimensional hydrodynamic solver;

[0013] Step 4: Establish a dry-wet dynamic boundary processing module with memory and compensation mode, make it work synchronously with the hydrodynamic solver, and together form a complete two-dimensional hydrodynamic model, which has the ability to process the dynamic changes of water areas in the flood propagation of braided complex rivers in real time;

[0014] Step 5: For the target section of the braided river, set its inlet and outlet flow and water level open boundary conditions, carry out simulation of the plane flow field of the braided river, save the water level distribution and flow field results calculated by the model to the hard disk, and display them on the computer screen to provide feedback to the user.

[0015] In step 1, the two-dimensional computational grid of the river is prepared as follows:

[0016] Step 11: Referring to the waterside line and beach edge line in the measured topographic map of the braided river channel, the river channel is divided into three types of areas. The first is the beach in the river channel, which is the area between the beach edge line of the beach on the left and right inner sides of the river channel and the outer boundary of the river channel, marked with -1; the second is the sandbar in the river channel, which is submerged only when the flow rate is large. The sandbar in the river channel is marked with non-repeating positive integers such as 1, 2, 3, ..., N, where N is the total number of sandbars; the third is the channel between the beach and the sandbar, marked with 0;

[0017] Step 12: Draw the dike line according to the selected spatial step length, and take the area within the dike as the calculation area; draw the beach-channel boundary line along the beach edge line or waterside line in the river channel according to the selected spatial step length;

[0018] Step 13: Generate the surface areas of the side beach, groove flow channel, and center beach respectively, and obtain the computational grids of each partition of the three types of areas by generating each surface area separately.

[0019] In step 2, the terrain of the interpolation grid point is specifically:

[0020] Step 21: Corresponding to the classification of the calculation area in step 11, the measured terrain scattered points located in the beach, ditch, and sandbar divisions are marked with -1, 0, and positive integers respectively;

[0021] Step 22: Pair the grid points with the terrain scattered points through identification to form a "grid point-terrain scattered point" grouping;

[0022] Step 23: Within each group, the terrain of the grid points is interpolated based on the elevation data of the scattered terrain points.

[0023] The method of two-dimensional computational grid division of the braided river channel in step 1 and group terrain interpolation of the grid of the braided river channel in step 2 is referred to as the "partition division and group interpolation" fine modeling method. Its advantage is that the river channel and riverbed topography established based on this method can smoothly fit the flow path boundary, and the topographic measurement points above and below the water will not interfere with each other when interpolating the grid point terrain, thereby achieving a two-dimensional modeling effect of the river channel with smooth and simple groove boundaries, clear and accurate flow paths and river conditions, and accurate depiction of the braided river channel morphology, effectively avoiding the occurrence of groove boundary fragmentation, flow path blockage and other phenomena.

[0024] In step 3, the construction of the basic two-dimensional hydrodynamic solver is specifically as follows:

[0025] Step 31: Use the depth-averaged two-dimensional shallow water equations as the governing equations;

[0026] Step 32: Use the θ semi-implicit method to discretize the water level gradient term in the governing equations to eliminate the limitation on the computational stability of the solver caused by fast surface gravity waves;

[0027] Step 33: Use the Euler-Lagrange method to solve the convective terms in the governing equations so that the time step of the solver's calculation is not restricted by the stability conditions related to the grid scale;

[0028] Step 34: Use the finite volume method to discretize the governing equations to ensure the conservation of water volume in hydrodynamic calculations;

[0029] Step 35: Realize the coupled calculations of the water level gradient and the flow field in the river channel area by constructing and iteratively solving an algebraic system of equations based on velocity-pressure coupling.

[0030] The hydrodynamic solver belongs to the type and category of implicit two-dimensional hydrodynamic models. This type of solver can be applied to the solution of various subcritical, supercritical, and complex water flows and has the advantages of being able to use large time steps and high efficiency, and can solve the problem of instability in two-dimensional hydrodynamic calculations of braided channels. The main calculation parameters of the hydrodynamic solver include the channel roughness (n m ), etc. The operation mode of the hydrodynamic solver is: select a time step Δt and deduce the evolution process of the water flow step by step.

[0031] In step 4, the establishment of a dry-wet moving boundary processing module with a memory and compensation method, which coordinates and works synchronously with the hydrodynamic solver and jointly forms a complete two-dimensional hydrodynamic model, is specifically as follows:

[0032] Step 41: Detect the "wet → dry" state transition events and negative water depth cells of the unit, and update the dry-wet identification of the unit. Traverse all the cells that are wet at the beginning of this time step in the computational grid in ascending order of cell number. According to the latest cell water level solved by the hydrodynamic solver in this time step, calculate the water depth at the center of the cell, denoted by h. If h ≥ h0, the cell continues to remain in the wet state; if h < h0, it indicates that the cell has undergone a "wet → dry" state transition in this time step, and mark it as dry. For all the cells that have undergone a "wet → dry" transition in this time step, make a further judgment based on the cell water depth: if h > 0, the cell is an ordinary dry cell, and directly jump to step 43; if h < 0, mark the cell as a negative water depth cell, and mark the adjacent wet cells as affected cells, and execute step 42;

[0033] Step 42: Calculate the water overdraft and water return of the newly generated negative water depth dry unit, traverse all the newly generated negative water depth dry units according to the unit number, and execute steps 421 to 423;

[0034] Step 421: for a newly generated negative water depth dry unit, use its plane area A and water depth h to calculate its overdrawn water volume Ah at this time step;

[0035] Step 422: sort out and define the priorities of the affected units around the dry unit with negative water depth; at the same time, calculate and record the amount of water M that each of the affected units needs to return according to the water volume weights;

[0036] Step 423: specify the amount of water ETAC that the affected unit needs to return at each time step, and divide M by ETAC to obtain the number of water return times CI for each affected unit;

[0037] Step 43: The affected unit provides water flow according to the specified time step quota ETAC. For a certain affected unit, its water level decreases by ETAC in each time step. While the water level correction is being performed, its CI is reduced by 1. Once it is found that the CI has been reduced to 0, the "affected unit" identity of the unit is immediately eliminated, and all affected units are traversed according to the unit number to perform the above unit water level correction operation;

[0038] Step 44: Detect the unit "dry→wet" conversion event and update the unit identification. Traverse all dry units in the order of unit numbers to determine whether they will become wet due to the hydraulic conditions of neighboring units. For ordinary dry units, when the unit recovery flow condition is met, set its dry-wet flag to wet and the water depth to h0, and let it participate in the model calculation of the next time step. For negative water depth dry units, when the unit recovery flow condition is met, add a water level increment Δh to the memory variable η of the unit. When it is found that the updated η is higher than the bed surface by h0, immediately eliminate its "negative water depth unit" identity and set its dry-wet flag to wet. After this detection and update, all remaining dry units do not participate in the model calculation of the next time step;

[0039] Step 45: Based on the latest unit wetness information obtained in steps 41 to 44, update the wetness status of other grid elements such as nodes and edges in the calculation grid.

[0040] The wet-dry moving boundary processing module of the memory and compensation method described above needs to define the following variables before application. A wet-dry mark (variable DW) is defined for each cell of the grid to distinguish between wet cells and dry cells; at the same time, a negative water depth mark (variable NEG) is defined for each cell to distinguish between negative water depth dry cells and ordinary dry cells. The wet cells around the negative water depth cells are defined as "affected cells". According to the flow direction of each side of the negative water depth cell (refer to the state of the cell before it is converted from a wet cell to a negative water depth cell), the affected cells are divided into two categories: if the flow direction deviates from the negative water depth cell, it is called a receiving cell, otherwise it is called an ordinary neighboring cell. It is stipulated that the priority of the water return volume of the receiving cell is higher than that of the ordinary neighboring cell. For the negative water depth cell, its instantaneous water level (η) is used as a memory variable (Note: the negative water depth cell η is lower than the riverbed).

[0041] The amount of water returned by the affected unit at each time step is ETAC, which specifically refers to the height by which the water level of the water column in the two-dimensional unit is lowered due to the water expenditure correction at each time step.

[0042] The memory and compensation method for handling wet and dry dynamic boundaries has added a mechanism for handling negative water depth units compared to the traditional critical water depth method: the neighboring units of the negative water depth unit gradually return the overdrawn water volume of the negative water depth unit through step-by-step local correction, thereby maintaining the conservation of water volume calculated in the complex water flow model while ensuring stability. The new wet and dry dynamic boundary processing method solves the contradiction between the large time step hydrodynamic calculation and the small time step required to simulate the rapidly changing wet and dry boundary dynamics. The execution of the new wet and dry dynamic boundary processing method is independent of the hydrodynamic calculation and will not adversely affect the hydrodynamic calculation. In addition, this method is not only applicable to braided rivers (the most complex), but also to other types of rivers (such as straight, curved, and bifurcated rivers), and is a general simulation method.

[0043] Compared with the prior art, the present invention has the following beneficial effects:

[0044] The fine modeling through partitioning and group interpolation lays the foundation for accurately simulating the complex water flow plane flow field of braided rivers; the dry-wet dynamic boundary processing method using memory and compensation methods ensures the conservation of hydrodynamic calculations without affecting the stability, time step, execution process, and efficiency of hydrodynamic calculations; and finally solves the problems of large errors, poor effects, or low versatility of braided river water flow plane flow field simulation methods in the field of river simulation technology. The method of the present invention also has the advantages of high reliability, high calculation efficiency, and strong versatility. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 It is a schematic diagram of the method flow for implementing the present invention.

[0046] Figure 2It is a schematic diagram of the execution flow of the dry-wet dynamic boundary processing method implementing the memory and compensation mode of the present invention.

[0047] Figure 3 It is a schematic diagram of the classification and identification of planar areas of braided river channels according to the present invention.

[0048] Figure 4 It is a schematic diagram of the partitioning of the two-dimensional computational grid of a braided river channel according to the present invention.

[0049] Figure 5 It is a schematic diagram of the result of terrain interpolation of braided river channel calculation grid grouping according to the present invention.

[0050] Figure 6 It is a schematic diagram of the water boundary and flow field of the Niyang River in accordance with the present invention. DETAILED DESCRIPTION

[0051] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0052] As attached Figure 1 As shown, a method for simulating the plane flow field of braided river flow includes the following steps:

[0053] Step 1: Create a two-dimensional computational grid for the river channel, including classifying and marking the plane area of ​​the braided river channel, drawing the boundary line inside and outside the river channel, and generating computational grids by partition;

[0054] Step 2: interpolating the terrain of the grid points, including classifying and marking the scattered points and grid points of the braided river terrain, pairing the grid points with the scattered points of the terrain, and interpolating the terrain of the grid points, thereby completing the fine modeling of the riverbed topography and geomorphology in the planar area of ​​the braided river;

[0055] Step 3: numerically discretize the two-dimensional hydrodynamic control equations on a two-dimensional computational grid and construct a basic two-dimensional hydrodynamic solver;

[0056] Step 4: Establish a dry-wet dynamic boundary processing module with memory and compensation mode, make it work synchronously with the hydrodynamic solver, and together form a complete two-dimensional hydrodynamic model, which has the ability to process the dynamic changes of water areas in the flood propagation of braided complex rivers in real time;

[0057] Step 5: For the target reach of the braided river channel, set the inlet and outlet flow rates and water level open boundary conditions, conduct the simulation of the planar flow field of the braided river channel, save the water level distribution and flow field results calculated by the model to the hard disk, and display them on the computer screen for feedback to the user.

[0058] Figure 2 It is a schematic diagram of the execution process of the wet-dry moving boundary treatment method for the memory and compensation method of implementing the present invention. Figure 2 In

[0059] Step 41: Detect the "wet→dry" state conversion event and negative water depth units of the detection unit, and update the wet-dry identification (variable) of the unit. Traverse all the units that are wet at the beginning of this time step in the computational grid in ascending order of unit number, and calculate the water depth at the center of the unit (denoted by h) according to the latest unit water level solved by the hydrodynamic solver in this time step. If h≥h0, the unit continues to remain in the wet state; if h<h0, it indicates that the unit has undergone a "wet→dry" state transition in this time step, and mark it as dry. For all units that have undergone a "wet→dry" conversion in this time step, further judgment is made according to the unit water depth: if h>0, the unit is an ordinary dry unit, and directly jump to Step 43; if h<0, mark the unit as a negative water depth unit, mark the adjacent wet units as affected units, and execute Step 42.

[0060] Step 42: Carry out the calculation of the water volume overdraw and the returned water volume of the newly generated negative water depth dry units. All the newly generated negative water depth dry units can be traversed according to the unit number, and the following sub-steps are executed: ① For a newly generated negative water depth dry unit, use its planar area A and water depth h to calculate the water volume overdrawn by it in this time step Ah; ② Sort out and specify the priorities of the affected units around the negative water depth dry units, and calculate and record the water volumes M that they need to return respectively according to the water volume weights of the affected units; ③ Specify the water volume (ETAC) that the affected units need to return in each time step, and divide M by ETAC to obtain the water volume return times (CI) of each affected unit.

[0061] Step 43: The affected units give out water according to the specified time step quota (ETAC). For example, for a certain affected unit, its water level drops by ETAC in each time step; while performing the water level correction, its CI decreases by 1. Once it is found that CI has decreased to 0, immediately eliminate the "affected unit" status of the unit. All the affected units can be traversed according to the unit number to perform the above unit water level correction operation.

[0062] Step 44: Detect the unit "dry→wet" conversion event and update the unit identification. Traverse all dry units in the order of unit numbers to determine whether they will become wet due to the hydraulic conditions of neighboring units. For ordinary dry units, when the unit recovery flow conditions are met, set its dry-wet flag to wet and the water depth to h0, and let it participate in the model calculation of the next time step. For negative water depth dry units, when the unit recovery flow conditions are met, add a water level increment Δh to the memory variable (η) of the unit. When it is found that the updated η is higher than the bed surface by h0, its "negative water depth unit" identity is immediately eliminated, and its dry-wet flag is set to wet. After this detection and update, all remaining dry units do not participate in the model calculation of the next time step.

[0063] Step 45: Based on the latest unit wetness information obtained in steps 41 to 44, update the wetness status of other grid elements such as nodes and edges in the calculation grid.

[0064] Figure 3 It is a schematic diagram of the classification and identification of planar areas of braided rivers according to the present invention (taking the estuary section of the Niyang River, a first-level tributary of the Yarlung Zangbo River, as an example). Figure 3 middle,

[0065] According to the waterside line and beach edge line in the measured topographic map (1:10000), the target river section is divided into three types of areas: ① Beach, the area between the beach edge line of the beach on the left and right inner sides of the river and the outer boundary of the river (including beach, highland, etc.), which has a large elevation and does not flow under small and medium flow, and is marked with -1; ② Sandbars in the river, which are only submerged when the flow is large, are marked with non-repeating numbers, such as 1, 2, 3, ..., N (N is the total number of sandbars); ③ The channel flow between the beach and the sandbar is marked with 0. Correspondingly, the measured terrain scattered points located in the beach, ditch, and sandbar divisions are marked with -1, 0, and positive integers respectively. On this basis, the grid points are paired with the terrain scattered points through marking to form several "grid point-terrain scattered point" groups.

[0066] Figure 4 It is a schematic diagram of the two-dimensional computational grid partitioning of a braided river channel according to the present invention (taking the estuary section of the Niyang River, a first-level tributary of the Yarlung Zangbo River, as an example). Figure 4 middle,

[0067] According to the set grid size (50m in this example), Figure 3 The three types of regions in the grid are meshed separately and appropriately encrypted in the trench range. The resulting triangular unstructured grid contains 91396 cells and 46433 nodes. Then, the cell center is selected as the grid point of "control terrain", and the cell center in each partition is marked with the corresponding partition mark to prepare for group terrain interpolation.

[0068] Figure 5 It is a schematic diagram of the result of the grid grouping terrain interpolation of the braided river channel calculation according to the present invention (taking the estuary section of the Niyang River, a first-level tributary of the Yarlung Zangbo River, as an example). Figure 5 middle,

[0069] As can be seen from the figure, when using the conventional grid division and terrain interpolation method (using uniform grids, mixing the above-water and underwater terrain scatter points for grid point terrain interpolation), the two-dimensional grid cannot smoothly fit the flow path boundary, and the above-water and underwater terrain measurement points interfere with each other when interpolating the grid point terrain, the resulting trench water boundary is very fragmented and it is difficult to clearly reflect the flow path of the braided river. When the partition division and group interpolation method is used (in each "grid point-terrain scatter point" group, the terrain of the grid point is interpolated based on the elevation data of the terrain scatter points), the resulting trench boundary is smooth and simple, which can clearly depict the flow path and river flow of the braided river, and the modeling effect is consistent with the real river morphology.

[0070] Figure 6 It is a schematic diagram of the water boundary and flow field of the Niyang River calculated by implementing the method of the present invention (taking the estuary section of the Niyang River, a first-level tributary of the Yarlung Zangbo River, as an example). Figure 6 middle,

[0071] By comparison, it is found that the dry season flow path ( Figure 6 b) Compared with the measured dry season flow path ( Figure 3 ) are in good agreement, while the flow area calculated by the two-dimensional model using the traditional critical water depth method is much larger than the measured dry season water area (i.e., the calculated plane flow field is seriously distorted). In the upper section of the Niyang River in the calculation area, the water area and flow field simulated by the two methods are relatively close. However, the flow field calculated by the two-dimensional model using the traditional critical water depth method shows a significant enhancement after the water flow passes through the area with complex river flow units, and when this water flow enhancement phenomenon develops to the downstream estuary area, many grooves that are not actually flowing also become flow channels. These scenarios are significantly different from the measured flow paths.

[0072] The water conservation error calculated by the two-dimensional hydrodynamic model shows that the flow rate at the outlet of the study area calculated by the model using the traditional critical water depth method is almost doubled relative to the total inflow flow, and the water conservation error is 94.31%; in contrast, the water conservation error of the model using the method of the present invention when simulating braided rivers is less than 1.0%, ensuring the accuracy of the two-dimensional hydrodynamic calculation results.

[0073] The above embodiments are merely examples of the technical solutions of the present invention. The inversion model construction method involved in the present invention is not limited to the contents described in the above embodiments, but is subject to the scope defined in the claims. Any modification, supplement or equivalent replacement made by a person skilled in the art of the present invention on the basis of the embodiment is within the scope of protection claimed by the present invention.

Claims

1. A method for simulating the planar flow field of braided river flow, characterized in that: The steps include: Step 1: Create a two-dimensional computational grid for the river channel, including classifying and marking the plane area of ​​the braided river channel, drawing the boundary line inside and outside the river channel, and generating computational grids by partition; Step 2: interpolating the terrain of the grid points, including classifying and marking the scattered points and grid points of the braided river terrain, pairing the grid points with the scattered points of the terrain, and interpolating the terrain of the grid points, thereby completing the fine modeling of the riverbed topography and geomorphology in the planar area of ​​the braided river; Step 3: numerically discretize the two-dimensional hydrodynamic control equations on a two-dimensional computational grid and construct a basic two-dimensional hydrodynamic solver; Step 4: Establish a dry-wet dynamic boundary processing module with memory and compensation mode, make it work synchronously with the hydrodynamic solver, and together form a complete two-dimensional hydrodynamic model, which has the ability to process the dynamic changes of water areas in the flood propagation of braided complex rivers in real time; Step 5: For the target section of the braided river, set its inlet and outlet flow and water level open boundary conditions, carry out simulation of the plane flow field of the braided river, save the water level distribution and flow field results calculated by the model to the hard disk, and display them on the computer screen to provide feedback to the user.

2. A method for simulating the planar flow field of braided river flow according to claim 1, characterized in that: In step 1, the two-dimensional computational grid of the river is prepared as follows: Step 11: Referring to the waterside line and beach edge line in the measured topographic map of the braided river channel, the river channel is divided into three types of areas. The first is the beach in the river channel, which is the area between the beach edge line of the beach on the left and right inner sides of the river channel and the outer boundary of the river channel, marked with -1; the second is the sandbar in the river channel, which is submerged only when the flow rate is large. The sandbar in the river channel is marked with non-repeating positive integers such as 1, 2, 3, ..., N, where N is the total number of sandbars; the third is the channel between the beach and the sandbar, marked with 0; Step 12: Draw the dike line according to the selected spatial step length, and take the area within the dike as the calculation area; draw the beach-channel boundary line along the beach edge line or waterside line in the river channel according to the selected spatial step length; Step 13: Generate the surface areas of the side beach, groove flow channel, and center beach respectively, and obtain the computational grids of each partition of the three types of areas by generating each surface area separately.

3. A method for simulating the planar flow field of braided river flow according to claim 2, characterized in that: In step 2, the terrain of the interpolation grid point is specifically: Step 21: Corresponding to the classification of the calculation area in step 11, the measured terrain scattered points located in the beach, ditch, and sandbar divisions are marked with -1, 0, and positive integers respectively; Step 22: Pair the grid points with the terrain scattered points through identification to form a "grid point-terrain scattered point" grouping; Step 23: Within each group, the terrain of the grid points is interpolated based on the elevation data of the scattered terrain points.

4. A method for simulating the planar flow field of braided river flow according to claim 1, characterized in that: In step 3, the basic two-dimensional hydrodynamic solver is constructed as follows: Step 31: Use the depth-averaged two-dimensional shallow water equation as the control equation; Step 32: Use the θ semi-implicit method to discretize the water level gradient term in the control equation to eliminate the limitation on the solver calculation stability caused by the fast surface gravity wave; Step 33: Use the Euler-Lagrangian method to solve the convection term in the control equation so that the time step of the solver calculation is not limited by the stability condition related to the grid size; Step 34: Discretize the governing equations using the finite volume method to ensure the conservation of water volume in hydrodynamic calculations; Step 35: Achieve the coupled calculations of water level gradients and flow fields in the river channel area by constructing and iteratively solving an algebraic equation system based on velocity-pressure coupling.

5. A method for simulating the planar flow field of braided river flow according to claim 1, characterized in that: In the said Step 4, a dry-wet moving boundary processing module with a memory and compensation method is established to work coordinately and synchronously with the hydrodynamic solver, and they are jointly formed into a complete two-dimensional hydrodynamic model, specifically: Step 41: Detect the "wet → dry" state transition events and negative water depth cells, and update the wet-dry identifiers of the cells. Traverse all the cells that are wet at the beginning of this time step in the computational grid in ascending order of cell numbers. According to the latest cell water level solved by the hydrodynamic solver at this time step, calculate the water depth at the cell center, denoted as h. If h ≥ h0, the cell continues to remain in the wet state; if h < h0, it indicates that the cell has undergone a "wet → dry" state transition at this time step, and mark it as dry. For all the cells that have undergone a "wet → dry" transition at this time step, make a further judgment based on the cell water depth: if h > 0, this cell is an ordinary dry cell, and directly jump to Step 43; if h < 0, mark this cell as a negative water depth cell, and mark the adjacent wet cells as affected cells, and execute Step 42; Step 42: Conduct the calculation of water volume overdraw and water return volume for the newly generated negative water depth dry cells. Traverse all the newly generated negative water depth dry cells in ascending order of cell numbers, and execute Steps 421 to 423; Step 421: For a newly generated negative water depth dry cell, use its planar area A and water depth h to calculate the overdrawn water volume Ah at this time step; Step 422: Sort out and specify the priorities of the affected cells around the negative water depth dry cells; at the same time, calculate and record the respective water volumes M that they need to return according to the water volume weights of the affected cells; Step 423: Specify the water volume ETAC that the affected cells need to return at each time step, and divide M by ETAC to obtain the water return times CI of each affected cell; Step 43: The affected cells give out water according to the specified time step quota ETAC. For a certain affected cell, its water level drops by ETAC at each time step; while performing the water level correction, its CI decreases by 1. Once it is found that CI has decreased to 0, immediately eliminate the "affected cell" status of this cell. Traverse all the affected cells in ascending order of cell numbers and perform the above-mentioned cell water level correction operation; Step 44: Detect the unit "dry→wet" conversion event and update the unit identification. Traverse all dry units in the order of unit numbers to determine whether they will become wet due to the hydraulic conditions of neighboring units. For ordinary dry units, when the unit recovery flow condition is met, set its dry-wet flag to wet and the water depth to h0, and let it participate in the model calculation of the next time step. For negative water depth dry units, when the unit recovery flow condition is met, add a water level increment Δh to the memory variable η of the unit. When it is found that the updated η is higher than the bed surface by h0, immediately eliminate its "negative water depth unit" identity and set its dry-wet flag to wet. After this detection and update, all remaining dry units do not participate in the model calculation of the next time step; Step 45: Based on the latest unit wetness information obtained in steps 41 to 44, update the wetness status of nodes, edges and other grid elements in the computational grid.

Images