A method and system for constructing explicit unstructured two-dimensional hydrodynamic models

Abstract

This invention relates to a method and system for constructing an explicit unstructured two-dimensional hydrodynamic model. The method includes constructing an unstructured mesh and identifying wet and dry states based on critical water depths, adapting to complex boundaries; discretizing the governing equations using the finite volume method based on conserved two-dimensional shallow water equations, and simultaneously solving the pressure and velocity terms through an explicit pressure-velocity coupling mechanism to achieve flow field prediction and correction and water depth updates; and combining optimization techniques such as adaptive time stepping and vectorization acceleration to perform time-iterative calculations and output visualized results. Its advantages include: effectively suppressing numerical oscillations under explicit in-situ meshes without requiring implicit or iterative algorithms, thus improving simulation accuracy and stability; improving computational efficiency through algorithm optimization, supporting large-scale mesh simulations; and adapting to complex boundaries such as estuaries and coastlines, providing a scientific basis for integrated estuarine and coastal management and intertidal ecological protection and restoration.

Description

Technical Field

[0001] This invention belongs to the field of water conservancy engineering and relates to the field of hydrodynamic numerical simulation technology, especially the method and system for constructing explicit unstructured two-dimensional hydrodynamic models. Background Technology

[0002] Hydrodynamic numerical simulation is a core technology for scientific decision-making in hydraulic and ocean engineering, playing a crucial role in port engineering, coastal engineering, flood control and disaster reduction, and watershed management. With the increasing sophistication of engineering requirements, unstructured meshes have gradually become the preferred solution for high-precision hydrodynamic simulation due to their adaptability to complex boundaries (such as estuaries, intertidal zones, and islands).

[0003] Current mainstream hydrodynamic numerical simulations are mainly based on implicit coupling. By constructing large sets of equations to solve pressure and velocity terms simultaneously, numerical stability can be guaranteed. However, this requires large-scale matrix inversion operations, resulting in extremely low computational efficiency and making it difficult to meet the needs of real-time simulation of large-scale watersheds.

[0004] On the other hand, although explicit coupling methods in traditional numerical simulations avoid matrix operations, they have significant drawbacks in unstructured grids: the decoupling of pressure and velocity fields can easily lead to numerical oscillations, especially in high velocity gradient regions (such as sharp bends in rivers and scour zones below dams), where simulation accuracy drops sharply; in order to maintain stability, extremely small time steps are required, which offsets the efficiency advantages of explicit algorithms.

[0005] Furthermore, existing two-dimensional hydrodynamic models are mostly based on structured grid designs, which are difficult to adapt to the irregular cell shapes of unstructured grids, and boundary interpolation errors further amplify the simulation deviation.

[0006] Therefore, how to construct a two-dimensional hydrodynamic model in an unstructured grid system that combines high accuracy, high stability, and high efficiency, and can realize an explicit pressure-velocity coupling mechanism, has become a key issue in overcoming the application limitations of unstructured two-dimensional hydrodynamic models. Summary of the Invention

[0007] The purpose of this invention is to solve the above-mentioned problems existing in the prior art and to provide a method and system for constructing an explicit unstructured two-dimensional hydrodynamic model.

[0008] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0009] A method for constructing an explicit unstructured two-dimensional hydrodynamic model includes the following steps:

[0010] S1: Mesh construction and data preparation. Construct an unstructured mesh under ideal working conditions, identify the dry and wet states of the mesh based on the critical water depth threshold and dynamically process the interface flux, and adapt complex boundary conditions according to the flux characteristics of the mesh boundary edges.

[0011] S2: A two-dimensional shallow water dynamic model is constructed based on the conservation type two-dimensional shallow water equation. The control equation is discretized using the finite volume method. The pressure and velocity terms in the continuity equation and momentum equation are solved simultaneously through an explicit pressure-velocity coupling mechanism to achieve flow field prediction, correction and water depth update.

[0012] S3: Set the model boundary conditions and initial conditions. The boundary conditions include Dirichlet boundary, Neumann boundary or absorbing boundary, and the initial conditions include topography, water depth and flow velocity parameters.

[0013] S4: An explicit time-stepping method is used for time iteration calculation. The time step size is dynamically adjusted based on the flow field characteristics, and the iteration continues until the set total number of time steps is reached to ensure numerical stability.

[0014] S5: Outputs the hydrodynamic field simulation results at each time point, visualizes the flow rate and velocity parameters, and analyzes the characteristics of water flow.

[0015] Preferably, the identification of the wet / dry state of the grid in step S1 includes: introducing a critical water depth threshold as a judgment variable; when the water depth within the grid is ≥ the critical threshold, it is determined to be a wet grid; otherwise, it is determined to be a dry grid. Interface flux processing includes the following scenarios:

[0016] (1) When both sides of the grid are dry grids, or when the liquid level of the wet grid is lower than that of the dry grid, there is no flux at the interface;

[0017] (2) When both sides of the grid are wet grids, or when the liquid level of the wet grid is greater than or equal to that of the dry grid, there is flux at the interface.

[0018] Preferably, the adaptation to complex boundary conditions in step S1 includes: when all three edges of the grid are non-flux boundaries, the grid does not participate in the calculation; when all three edges of the grid are flux boundaries, the grid is calculated normally; when some of the three edges of the grid are flux boundaries, the grid only calculates the mass equation and does not calculate the momentum equation.

[0019] Preferably, the finite volume method discretization in step S2 includes: transforming the volume integral of the governing equation into the boundary integral using the Gaussian divergence theorem, and obtaining the flow rate update formula after discretization.

[0020]

[0021] Water depth update formula

[0022]

[0023] Among them, Q i A is the flow vector of unit i; i F is the area of ​​cell i; conv F pres Fdiff These represent the convection flux, pressure flux, and diffusion flux, respectively; Δt is the time step; h i Let Q be the water depth of unit i. f Let n be the flow of the boundary edge. f Let L be the edge normal vector. f Let be the side length.

[0024] Preferably, the explicit pressure-velocity coupling mechanism in step S2 employs an explicit coupling algorithm based on an improved Rhie-Chow interpolation. A pressure gradient-related dissipation term is added to the unit interface velocity obtained through linear interpolation to correct the flow field. The equation for correcting the flow field is as follows:

[0025]

[0026] in, All are obtained through linear interpolation, as shown in the following formula:

[0027]

[0028] e CF Let be the unit direction vector of CF.

[0029] Preferably, in step S2, the convection term in the momentum equation is treated using the upwind mode. The direction of the water flow is determined by calculating the dot product of the edge velocity and the edge unit normal vector. The formula for calculating the edge velocity is as follows:

[0030]

[0031] Among them, Q C and Q Fi Let the flow rates be at unit center C and adjacent unit F, respectively. The judgment equation is as follows:

[0032]

[0033] ε is the direction determination vector. When ε > 0, the water flow direction is from the center cell C to the adjacent cell Fi; when ε < 0, the water flow direction is from the adjacent cell Fi to the center cell C. Among these, Q in the convection term... C The following formula will be used for assignment:

[0034]

[0035] The relevant equations for flux calculation are as follows:

[0036]

[0037] F conv F pres F diffThese are the convection flux, pressure flux, and diffusion flux, respectively, h fi ν is the water depth at the interface. t L is the viscosity coefficient. f Let A be the length of the common edge between the interface of cell center C and adjacent cell F. i Let C be the area of ​​the unit center.

[0038] Preferably, the dynamic adjustment of the time step in step S4 is based on the CFL condition, and the maximum allowable time step is dynamically calculated by combining the local wave velocity and the minimum characteristic length of the grid cell. When the velocity gradient is large or the grid cell is small, the time step is automatically reduced, and vice versa.

[0039] Preferably, in step S2, NumPy vectorization operations and Numba library just-in-time compilation technology are used to optimize the core calculation process. The nested loop logic of cell boundary value calculation and flux accumulation is rewritten into matrix operations, and the @numba.jit decorator is added to the core functions of gradient reconstruction and time step iteration for machine code compilation.

[0040] Preferably, step S1 also includes an adaptive mesh refinement function, which identifies high-variable regions in the flow field based on the velocity gradient modulus or the water depth gradient modulus, and further subdivides mesh cells whose gradient values ​​exceed a set threshold, dividing the triangular cells into 4 sub-triangular cells.

[0041] Another objective of this invention is to provide an explicit unstructured two-dimensional hydrodynamic model construction system.

[0042] To achieve the second objective mentioned above, the technical solution adopted by the present invention is as follows:

[0043] An explicit unstructured two-dimensional hydrodynamic model construction system, comprising:

[0044] The mesh and data processing module is used to construct unstructured meshes, identify the wet and dry states of meshes based on critical water depth thresholds, and adapt complex boundaries according to the flux characteristics of boundary edges.

[0045] The model building module is used to predict, correct and update the flow field and water depth based on the conserved two-dimensional shallow water equations, through the finite volume method to discretize the control equations and the explicit pressure-velocity coupling mechanism.

[0046] The parameter configuration module is used to set boundary conditions, initial conditions, and dynamic time step parameters.

[0047] The computation iteration module is used for dynamic iterative computation using an explicit time-stepping method;

[0048] The results output module is used to output simulation results and perform visualization processing and water flow motion characteristic analysis.

[0049] Due to the adoption of the above technical solution, the beneficial effects obtained by the present invention include:

[0050] 1. This invention proposes a method for solving pressure and velocity fields coupled based on explicit isotopic grids. It eliminates the need for implicit or iterative algorithms and adds a dissipative term to the interface velocity using an improved Rhie-Chow interpolation technique. This solves the numerical oscillation problem that traditional explicit methods are prone to in unstructured grids (especially in high velocity gradient regions such as sharp bends in rivers and scour zones downstream of dams). At the same time, it significantly improves the efficiency of coupled solution while ensuring computational accuracy, avoids large-scale matrix inversion operations, and reduces hardware requirements.

[0051] 2. This invention introduces variables to dynamically determine the wet and dry state of the grid (wet grid: water depth ≥ threshold; dry grid: water depth < threshold), and adaptively processes flux based on the wet and dry combination of the grids on both sides of the interface (e.g., "one dry and one wet with the liquid level of the wet grid higher than that of the dry grid"). This solves the simulation instability problem caused by frequent wet and dry transitions of the grid during high and low tides, thereby accurately simulating the water level change process in complex water areas such as the intertidal zone and estuary, and improving the model's adaptability to wet and dry boundaries.

[0052] 3. This invention breaks through the efficiency and stability bottlenecks of traditional unstructured grid hydrodynamic models by using three core technologies: explicit coupling mechanism, dry and wet treatment, and complex boundary adaptation. Combined with adaptive algorithms, vectorization acceleration, and multi-scenario functional expansion, it forms a set of numerical simulation methods that take into account high accuracy, high efficiency and flexibility, which can be directly applied to engineering practices such as estuary and coastal integrated management and intertidal ecological restoration. Attached Figure Description

[0053] Figure 1 This is a flowchart illustrating the structural process of the explicit unstructured two-dimensional hydrodynamic model construction method of the present invention.

[0054] Figure 2 This is a technical roadmap of an embodiment of the explicit unstructured two-dimensional hydrodynamic model construction method of the present invention.

[0055] Figure 3 This is a schematic diagram of a single mesh and its relationship with adjacent meshes in the unstructured two-dimensional hydrodynamic model of the present invention.

[0056] Figure 4 This is a comparison diagram of the effects of the pressure-velocity coupling algorithm before and after its application in an embodiment of the present invention.

[0057] Figure 5 This is a comparison chart of the simulation results of the model proposed in this embodiment and the Telemac2D model.

[0058] Figure 6 This is a flow velocity distribution diagram at a specific moment in an embodiment of the present invention. Detailed Implementation

[0059] Please see Figure 1-6 As shown, this invention mainly provides a method and system for constructing an explicit unstructured two-dimensional hydrodynamic model. Based on a Python-language unstructured two-dimensional hydrodynamic model, it establishes an efficient data interaction pathway between pressure and velocity fields through mesh topology mapping technology. By constructing an explicit coupling mechanism in the unstructured two-dimensional hydrodynamic model, it achieves the goal of efficiently and accurately simulating hydrodynamic processes in complex water bodies. During this process, by combining mature unstructured mesh discretization technology and a self-written explicit pressure-velocity coupling algorithm, Python and its scientific computing libraries are used to automatically complete model coupling, data transmission, and scientific computing, achieving the goal of high-precision and rapid simulation of two-dimensional water flow motion under complex boundary conditions. Based on this, the method specifically includes the following steps:

[0060] S1: Mesh construction and data preparation. Construct an unstructured mesh under ideal working conditions, identify the dry and wet states of the mesh based on the critical water depth threshold and dynamically process the interface flux, and adapt complex boundary conditions according to the flux characteristics of the mesh boundary edges.

[0061] In this embodiment, identifying the wet / dry state of the grid includes: introducing a critical water depth threshold as a judgment variable; when the water depth in the grid is ≥ the critical threshold, it is determined to be a wet grid; otherwise, it is determined to be a dry grid. The interface flux processing includes the following situations: (1) when both grids are dry grids, or when the liquid level of the wet grid is lower than that of the dry grid, there is no flux at the interface; (2) when both grids are wet grids, or when the liquid level of the wet grid is greater than or equal to that of the dry grid, there is flux at the interface.

[0062] In this embodiment, by using critical water depth variables and flux processing rules, the dry and wet states of the grid are accurately identified, solving the problem of dry-wet transition in scenarios such as high tide and low tide, avoiding interface flux calculation errors, and improving the model's simulation accuracy for shallow water areas (such as the intertidal zone). At the same time, the grid is processed according to different boundary flux characteristics, reducing interpolation errors caused by complex boundaries (such as shorelines and buildings), ensuring the conservation of momentum and mass in the boundary region, and improving the model's adaptability to irregular terrain.

[0063] Specifically, it is as follows:

[0064] 1. When generating the mesh, some meshes undergo wet-dry transitions during tidal cycles. The numerical model needs to accurately identify the mesh's wet-dry state and perform certain processing on the interface flux. The specific steps for interface flux processing include:

[0065] 1.1 Introducing the variable h min To determine the wet / dry state of a mesh, if h > h within that mesh... min The mesh is determined to be wet when h ≤ h min It is determined to be a dry mesh;

[0066] 1.2 Both meshes are dry, and there is no flux at this interface;

[0067] 1.3 Both meshes on both sides are wet, and the flux at this interface is calculated normally;

[0068] 1.4 When one mesh is dry and the other is wet, and the liquid level in the wet mesh is lower than that in the dry mesh, there is still no flux at the interface.

[0069] 1.5 When one grid is dry and the other is wet, and the liquid level in the wet grid is greater than or equal to that in the dry grid, there is flux at the interface.

[0070] 2. When constructing a computational grid adapted to complex boundaries, three cases need to be considered: all three edges of the grid have flux, none of the three edges have flux, and some edges of the grid have flux. The specific steps include:

[0071] 2.1 If all three edges of the grid are flux-free boundaries, the grid will not be included in the calculation.

[0072] 2.2 The three edges of the grid are all flux boundaries. The grid can be calculated normally and no special treatment is required.

[0073] 2.3. For grids where the three edges are partially flux-bearing, only the mass equation is calculated, not the momentum equation.

[0074] S2: A two-dimensional shallow water dynamic model is constructed based on the conservation type two-dimensional shallow water equation. The control equation is discretized using the finite volume method. The pressure and velocity terms in the continuity equation and momentum equation are solved simultaneously through an explicit pressure-velocity coupling mechanism to achieve flow field prediction, correction and water depth update.

[0075] The specific process is as follows:

[0076] 1. The unstructured two-dimensional hydrodynamic equations consider the relationship between flow rate (Q) and water depth (h), and their continuity equation is as follows:

[0077]

[0078] Where h is the water depth, t is the time, and Q is the flow vector (Q = hU). This represents the divergence operator.

[0079] 2. In the momentum equation of the unstructured two-dimensional hydrodynamic governing equations, to simplify the calculation, the source term in the formula can be ignored, and only the relationship between the time discrete term, convection term, and pressure term is considered. This simplification method reduces the amount of calculation while keeping the original monotonicity of the equation unchanged. The unstructured two-dimensional hydrodynamic governing equations under complex boundary conditions have a good fit. The equation is as follows:

[0080]

[0081] Where Q is the flow vector (Q = hU, where U is the velocity vector), and g is the gravitational acceleration. Let ν represent the tensor product, η be the water surface (η = h + z, where h is the water depth and z is the ground elevation), and ν be the water surface. t denoted as eddy viscosity coefficient, and u as flow velocity.

[0082] 3. The finite volume method is used to discretize the governing equations. The discretization process of the finite volume method revolves around "control volume integral - flux calculation - variable update," with the core being the transformation of the conserved two-dimensional shallow water equations into discrete algebraic equations for each grid cell. It uses the Gaussian divergence theorem to transform the "volume integral" into a "boundary integral," simplifying the discretization calculation. A schematic diagram is drawn, showing the relevant parameters for a single grid cell and its adjacent grid cells (as attached). Figure 3 As shown below:

[0083] In the diagram, C is the center cell, F1, F2, and F3 are the adjacent cells of cell C, f1, f2, and f3 are the common edges of the center cell and its adjacent cells, and h is the number of cells in the diagram. fi Q represents the water depth at the interface. fi n represents the flow at the interface. fi The boundary unit normal vector (the boundary of C and F2 is used as an example in the figure to show h) fi Q fi n fi ).

[0084] For any control volume Ω and its boundary Gauss's theorem is stated as:

[0085]

[0086] Where F is the flux vector.

[0087] Continuity equation (integral form):

[0088]

[0089] Where Ω represents the triangular cell region in the unstructured mesh. Let h be the boundary of the unit (three sides), t be the water depth, t be the time, Q be the flow vector, and n be the boundary unit normal vector.

[0090] Momentum equation (integral form):

[0091]

[0092] Where Q is the flow vector (Q = hU), Let g be the tensor product, g be the acceleration due to gravity, η be the water surface (η = h + z, where h is the water depth and z is the ground elevation), and ν be the tensor product.t denoted as eddy viscosity coefficient, and u as flow velocity.

[0093] Discretize the above integral equation into algebraic equations for each unit.

[0094] After discretization, the formula for updating the flow rate is:

[0095] Among them, Q i A is the flow vector of unit i; i F is the area of ​​cell i; conv F pres F diff These represent the convective flux, pressure gradient flux, and diffusion flux, respectively; Δt is the time step.

[0096] After discretization, the formula for updating the water depth can be obtained as follows:

[0097] Among them, h i Let A be the water depth of unit i. i Let Q be the area of ​​cell i. f Let n be the flow of the boundary edge. f Let L be the edge normal vector. f Let be the side length and Δt be the time step.

[0098] 4. When calculating the edge water depth, to facilitate the writing and understanding of the formulas, a schematic diagram is drawn showing the relevant parameters for a single grid and its adjacent grids (as shown in the attached diagram). Figure 2 (As shown). When calculating the water depth, first index the water depth of the center cell C and its adjacent cells Fi, and then use linear interpolation to calculate the water depth of the common edge fi of C and Fi, as shown in the following formula:

[0099] h fi =g Ci ·h C +g Fi ·h Fi

[0100]

[0101] Among them, h fi Let g be the water depth at the interface. Ci and g Fi d is the geometric interpolation factor corresponding to the adjacent cell Fi. Cf d Ff h represents the distance from the center cell C and its adjacent cell Fi to the intersection of the line connecting them and the common edge. C h Fi Let C be the water depth of the central cell C and its adjacent cell Fi.

[0102] 5. When calculating the water depth gradient, the Gauss-Green formula is used, and its formula is as follows:

[0103]

[0104] in, For the water level gradient, η fi S represents the water level at the interface. fi Let A be the boundary unit normal vector. i Let C be the area of ​​the unit center.

[0105] 6. Solve the momentum equation to predict the flow field.

[0106] The directionality of the convection process needs to be considered; therefore, an upwind approach is chosen to handle the convection process, and the dot product of the edge velocity and the edge unit normal vector is calculated. The edge velocity is calculated using linear interpolation, and the flow field is updated based on flux calculations. The equations for calculating the edge velocity are as follows:

[0107]

[0108] Among them, g c ,g F Q is a geometric factor. C and Q Fi Let the flow rates be at unit center C and adjacent unit F, respectively. The judgment equation is as follows:

[0109]

[0110] ε is the direction determination vector. The flow rate of the common edge of adjacent grid cells (see appendix) Figure 3 ), S fi Let Q be the normal vector at the interface between adjacent grid cells; when ε > 0, the water flow direction is from the central cell C to the adjacent cell Fi, and when ε < 0, the water flow direction is from the adjacent cell Fi to the central cell C; where Q is the normal vector at the interface between adjacent grid cells. C The following formula will be used for assignment:

[0111]

[0112] The relevant equations for flux calculation are as follows:

[0113]

[0114] F conv F pres F diff These are the convection flux, pressure flux, and diffusion flux, respectively, Q. c Q is the flow rate of the central cell C. Fi Q represents the flow rate of adjacent cell F. c Let Q be the flow vector of the central cell C.F S is the flow vector of adjacent grid F. fi h is the normal vector corresponding to the common edge of adjacent mesh cells. fi Let g be the water depth at the interface, g be the acceleration due to gravity, and η be the acceleration due to gravity. fi The water level at the interface. For the water level gradient between adjacent grid cells, ν t L is the viscosity coefficient. f Let A be the length of the common edge between the interface of cell center C and adjacent cell F. i Let C be the area of ​​the element center, and |CF| be the modulus of CF.

[0115] 7. In step S2, linear interpolation was used in processes such as water depth calculation. This interpolation method leads to decoupling of pressure and velocity at the element level, resulting in a checkerboard problem. We consider using an improved interpolation method based on Rhie-Chow interpolation. Rhie-Chow interpolation was originally used for pressure-velocity coupling in the Navier-Stokes equations and is generally used in implicit algorithms such as SIMPLE. It adds a dissipation term to the velocity obtained through linear interpolation on the element surface, equivalent to constructing a pseudo-momentum equation on the element surface. This simply mimics the small-scale pressure-velocity coupling pattern in a staggered mesh layout, closely linking the interface velocity to the pressure of adjacent elements, thus resolving the checkerboard problem. The relevant equations for correcting the flow field using the interpolation method in this paper are as follows:

[0116]

[0117] Among them, Q fi To correct the flow rate at the interface between adjacent grid cells, This represents the flow rate at the interface between adjacent grid cells after linear interpolation. To correct the pressure gradient between adjacent grid cells, This represents the pressure gradient between adjacent grid cells after linear interpolation. To correct for the flow field coefficients, the relevant parameters are calculated using the following formulas:

[0118]

[0119]

[0120] Among them, g C ,g Fi Let Q be the geometric factor, η be the flow rate, g be the gravitational acceleration, dt be the time step, and e be the geometric factor. CF Let |CF| be the unit direction vector of CF, and |CF| be the magnitude of CF.

[0121] S3: Combine different simulation scenarios to set model boundary conditions and initial conditions. Boundary conditions include Dirichlet boundary, Neumann boundary or absorbing boundary, and initial conditions include terrain, water depth and flow velocity parameters. Among them, combined with the characteristics of boundary elements, for boundary edges, the edge water depth is set to the water depth of the current element, the edge flow velocity is set to the flow velocity of the current element, and the edge flow rate is set to zero vector to reflect the constraint of boundary conditions on water flow.

[0122] S4: An explicit time-stepping method is used for time-iterative calculations. For each time step, functions are sequentially called to calculate the predicted flow field, correct the flow field, and update the water depth. The iteration process continues until the set total number of steps is reached, and the convergence of the calculation results is checked as needed to ensure numerical stability. The dynamic adjustment of the time step is based on the CFL condition, combined with the local wave velocity and the minimum characteristic length of the grid cells to dynamically calculate the maximum allowable time step. The time step is automatically reduced when the velocity gradient is large or the grid cells are small, and increased otherwise. The CFL number calculation formula is:

[0123]

[0124] Where i is the cell index within the computational domain, and j is the index of the neighboring cells surrounding the cell. R is the minimum feature length, which is the minimum distance from the center of all cells to the cell boundary.

[0125] S5: After the calculation is completed, output the hydrodynamic field simulation results at each time point, and visualize the simulation results in combination with the flow rate and velocity parameters to analyze the characteristics and laws of water flow.

[0126] In this embodiment, the least squares method is used instead of the existing summation method based on the edge vectors of a single cell in the gradient reconstruction stage of numerical computation. Specifically, for each grid cell, information from its 3-5 neighboring cells is integrated to construct a matrix based on the neighboring edge vectors and a function value difference vector. A weight matrix with the reciprocal of the cell spacing as the weight is introduced, and the gradient value is obtained by solving the least squares problem. This approach can effectively reduce interpolation errors under non-uniform grids, and is especially suitable for research areas with complex terrain and large differences in grid density, improving the stability and accuracy of gradient calculation.

[0127] In this embodiment, an adaptive time step mechanism is introduced into the time discretization scheme. Based on the CFL condition, the maximum allowable time step is dynamically calculated by combining the local wave velocity (or diffusion coefficient) and the minimum characteristic length of the grid cells, replacing the fixed time step setting. During the simulation, the local variation characteristics of the flow field are monitored in real time. When the velocity gradient is large or the grid cells are small, the time step is automatically reduced, and vice versa. This reduces unnecessary computational overhead and improves the overall simulation efficiency while ensuring numerical stability.

[0128] In this embodiment, NumPy vectorized operations and Numba library just-in-time compilation technology are used to optimize the core computation process. Nested loop logic such as cell boundary value calculation and flux accumulation is rewritten as matrix operations, leveraging NumPy's parallel computing capabilities to reduce Python loop overhead. Simultaneously, the `@numba.jit(nopython=True)` decorator is added to core functions such as gradient reconstruction and time step iteration for machine code compilation, significantly improving computation speed and enabling the model to efficiently handle large-scale simulation scenarios with more than 100,000 grid cells.

[0129] In this embodiment, the boundary condition processing module is extended to support flexible configuration of various boundary types. An additional boundary type file marks the properties of the mesh boundary edges (such as Dirichlet, Neumann, absorbing boundaries, etc.): for Dirichlet boundaries, the physical quantity values ​​at the boundary are directly given; for Neumann boundaries, a zero-gradient condition approximation is used; for absorbing boundaries, a sponge layer technique is introduced, setting an attenuation term in the boundary region to reduce wave reflection errors. Simultaneously, boundary parameters (such as given values, attenuation coefficients, etc.) are passed in through a configuration file, allowing adaptation to the boundary requirements of different simulation scenarios without modifying the core code.

[0130] In this embodiment, an adaptive mesh refinement function is added. High-variability regions in the flow field are identified based on gradient modulus (such as velocity gradient and water depth gradient). Mesh cells with gradient values ​​exceeding a set threshold are subdivided (e.g., a triangular cell is divided into four sub-triangles), and the mesh topology information and physical quantity interpolation are dynamically updated during the iteration process. By refining the mesh in high-precision regions and maintaining a sparse mesh in smooth regions, the simulation accuracy in key areas is ensured while controlling the overall computational load, balancing simulation efficiency and result reliability.

[0131] In this embodiment, the post-processing and verification module is enhanced. On the one hand, a comparison function with analytical solutions is added, which calculates quantitative indicators such as L2 error and maximum error between the simulation results and analytical solutions in real time, and displays the error distribution characteristics through charts to facilitate the evaluation of model accuracy. On the other hand, the visualization output format is expanded. In addition to generating animations, it supports exporting VTK format files, which can be imported into professional tools such as ParaView for advanced post-processing operations such as 3D flow field analysis and vector overlay. At the same time, visualization options for derivative physical quantities such as gradient vector field and flux distribution are added to provide richer data support for the analysis of water flow motion laws.

[0132] Based on the above, the present invention is demonstrated through the following examples:

[0133] like Figure 2 As shown, it specifically includes the following steps:

[0134] (1) Conduct benchmark tests and compare the simulation results before and after the application of the explicit pressure-velocity coupling algorithm.

[0135] (1-1) Construct a simple square computational domain, raise the local water level to allow water waves to propagate and reflect freely within the computational domain, and run the model code;

[0136] (1-2) Comparing the simulation results, it can be found that after applying the pressure-velocity coupling algorithm, the "chessboard oscillation" disappears significantly (e.g. Figure 4 Numerical stability is improved;

[0137] (2) Simulating the evolution of dam-break flow under complex terrain: mesh construction and data preparation

[0138] (2-1) Construct an unstructured mesh based on the computational domain and set the initial water depth and topography;

[0139] (2-2) Mark the mesh boundary type (entrance / exit / closed boundary);

[0140] (3) Input the parameters one by one into the two-dimensional shallow water dynamic model, perform flow field prediction and correction, and update the water depth.

[0141] (3-1) Dry and wet treatment: The model is enabled by default to automatically identify dry and wet grids. You can customize the critical water depth threshold (e.g., 0.01m) or determine the basis according to the simulation scenario (e.g., estuary, intertidal zone). When the water depth is lower than this value, the grid will automatically "dry up" and stop flux calculation.

[0142] (3-2) Finite Volume Method Discretization of Momentum Equations to Update Flow Field: The model uses input parameters and solves according to the set formula, discretizing the algebraic equations of each element. Update the grid flow field in the middle;

[0143] (3-3) Pressure-velocity coupled correction of flow field: The model automatically follows the formula Correcting the flow field eliminates the need for manual algorithm derivation;

[0144] (3-4) Finite volume method for discretizing continuity equations to update the water depth field: The model receives the corrected flow field data and discretizes it into the algebraic equations of each element according to the set formula. Update the grid depth field;

[0145] (3-5) Automatic Adaptive Calculation: Combining the model-corresponding algorithm, the evolution of the dam-break wave is simulated, and the dynamic change diagram of the flow field during the dam-break process is output (e.g., Figure 5 ).

[0146] (4) Results Output and Analysis

[0147] (4-1) After the model simulation is completed, output the water depth, flow velocity, flow rate and other data at each time point;

[0148] (4-2) Model visualization, and three-dimensional water level process plotting (e.g.) Figure 5 ), two-dimensional flow field plotting (such as Figure 6 );

[0149] (4-3) The simulation results are basically consistent with those of the mainstream open-source hydrodynamic model Telemac2D (e.g., Figure 5 The model employs a dam-break example proposed by Kawahara and Umetsu, which considers the evolution of the dam-break flow in initially dry but complex terrain. The model space is a rectangular region 75m long and 30m wide, and the terrain is defined as follows:

[0150]

[0151] Among them, z b The model uses a 50,000-grid grid, with an initial water depth of 1.875 m (x < 16 m). The remaining space is considered dry. The Manning coefficient is set to 0.018, and all boundaries use sliding wall boundary conditions. Compared with Telemac2D simulation results, the relative errors of global hydrodynamic parameters (water depth, velocity) are all less than 15%.

[0152] (4-4) Further analyze the flow field and explore the water flow process (e.g.) Figure 6 ).

[0153] Through the above steps, the unstructured two-dimensional hydrodynamic model based on the explicit pressure-velocity coupling algorithm can be run.

[0154] Another aspect of the present invention provides an explicit unstructured two-dimensional hydrodynamic model construction system, comprising:

[0155] (1) Mesh and data processing module, used to construct unstructured meshes, identify the dry and wet states of meshes based on critical water depth thresholds, and adapt complex boundaries according to the flux characteristics of boundary edges;

[0156] (2) Model building module, which is used to realize flow field prediction, correction and water depth update based on the conservation two-dimensional shallow water equation, through the finite volume method to discretize the control equation and the explicit pressure-velocity coupling mechanism;

[0157] (3) Parameter configuration module, used to set boundary conditions, initial conditions and dynamic time step parameters;

[0158] (4) Calculation iteration module, used for dynamic iterative calculation using explicit time stepping method;

[0159] (5) Result output module, used to output simulation results and perform visualization processing and water flow motion characteristic analysis.

[0160] In this embodiment, the modular design realizes an integrated process of mesh construction, model solving, parameter configuration and result output, supports flexible configuration of boundary conditions and calculation parameters, reduces the threshold for use, and facilitates algorithm optimization and function expansion.

[0161] It should be noted that this invention, by constructing an explicit unstructured two-dimensional hydrodynamic model, achieves efficient and stable simulation of water flow motion under complex boundary conditions while ensuring computational accuracy. Specific technical effects are as follows:

[0162] 1. An explicit pressure-velocity coupling mechanism is constructed based on the basic principles of fluid mechanics. The influence of unstructured mesh topology on the interaction between pressure and velocity fields is considered, which ensures the scientificity and accuracy of flow field simulation results under complex boundary conditions.

[0163] 2. The entire solution process is automated using Python and its scientific computing libraries. Calculation parameters and flow field boundary control are set through configuration files, which not only makes it easy for general technicians to get started quickly, but also allows for convenient algorithm optimization or functional expansion through modular code structure.

[0164] 3. The discrete solution of the hydrodynamic control equations is based on the mature finite volume method framework. The computational redundancy is reduced by the self-designed explicit coupling logic, and the mesh data processing and flow field visualization are realized by relying on third-party libraries. While maintaining a low computational time cost, the stability of the coupled solution of pressure and velocity fields under unstructured meshes is significantly improved, taking into account both simulation efficiency and result reliability.

[0165] In summary, this invention, through its technical system of "explicit coupling, efficient computation, stable discretization, and complex adaptation," balances simulation accuracy, computational efficiency, and adaptability to complex scenarios. It breaks through the limitations of traditional unstructured hydrodynamic models in terms of stability and engineering applications, providing a high-precision, low-cost numerical simulation tool for water conservancy and ocean engineering.

[0166] The foregoing descriptions and embodiments are provided to enable those skilled in the art to understand and apply the present invention. It will be apparent to those skilled in the art that various modifications can be easily made to these contents, and the general principles described herein can be applied to other embodiments without creative effort. Therefore, the present invention is not limited to the foregoing descriptions and embodiments. Improvements and modifications made by those skilled in the art based on the disclosure of the present invention without departing from its scope should be within the protection scope of the present invention.

Claims

1. A method for constructing an explicit unstructured two-dimensional hydrodynamic model, characterized in that, The method comprises the following steps: S1: mesh construction and data preparation, constructing an unstructured mesh under ideal working conditions, identifying the wet and dry states of the mesh based on the critical water depth threshold and dynamically processing the interface flux, and adapting the complex boundary conditions according to the flux characteristics of the mesh boundary edges; S2: constructing a two-dimensional shallow water dynamic model based on the two-dimensional shallow water equation, discretizing the control equation by using the finite volume method, synchronously solving the pressure term and the velocity term in the continuity equation and the momentum equation by using an explicit pressure-velocity coupling mechanism, realizing the prediction, correction and water depth update of the flow field, and the explicit pressure-velocity coupling mechanism adopts an explicit coupling algorithm improved based on the Rhie-Chow interpolation, and a dissipation term related to the pressure gradient is added to the unit interface velocity obtained by linear interpolation, so as to correct the flow field, and the equation for correcting the flow field is: where Q fi is the corrected flow rate at the interface of adjacent grid cells, is the flow rate at the interface of adjacent grid cells after linear interpolation, is the pressure gradient of adjacent grid cells after correction, is the pressure gradient of adjacent grid cells after linear interpolation, is the correction coefficient of the flow rate field, and the relevant parameter calculation formula is as follows: Among them, g C ,g Fi Let Q be the geometric factor, η be the flow rate, g be the gravitational acceleration, dt be the time step, and e be the geometric factor. CF Let |CF| be the unit direction vector of CF, and |CF| be the magnitude of CF. The convection term in the momentum equation is processed by using an upwind format, the dot product of the edge flow velocity and the unit normal vector of the edge is calculated to judge the flow direction, and the edge flow velocity calculation formula is: Among them, g c ,g F Q is a geometric factor. C and Q Fi Let the flow rates be at unit center C and adjacent unit F, respectively. The judgment equation is as follows: ε is a direction judgment vector, is the flow of the common edge of adjacent grid cells, S fi is the corresponding normal vector at the interface between adjacent grid cells; when ε > 0, the flow direction is that the center cell C flows to the adjacent cell F, and when ε < 0, the flow direction is that the adjacent cell F flows to the center cell C; wherein Q C will be assigned by the following formula: The flux calculation related equation is as follows: F conv , F pres , F diff are the convection, pressure and diffusion fluxes, respectively, Q c is the flow rate of the central cell C, Q Fi is the flow rate of the adjacent cell F, Q c is the flow rate vector of the central cell C, Q F is the flow rate vector of the adjacent cell F, S fi is the normal vector corresponding to the common edge of the adjacent cell, h fi is the water depth at the interface, g is the gravity acceleration, η fi is the water level at the interface, is the water level gradient of the adjacent cell, v t is the viscosity coefficient, L f is the length of the common edge of the interface between the cell center C and the adjacent cell F, A i is the area of the cell center C, |CF| is the length of CF; S3: setting the model boundary conditions and initial conditions, the boundary conditions include Dirichlet boundary, Neumann boundary or absorbing boundary, and the initial conditions include terrain, water depth and flow velocity parameters; S4: time iteration calculation is performed by using an explicit time stepping method, the time step is dynamically adjusted based on the flow field characteristics, and iteration is performed until a set total time step number is reached, so that the numerical stability is ensured; S5: outputting the water dynamic field simulation results at each time, visualizing the flow and flow velocity parameters and analyzing the flow motion characteristics.

2. The method of claim 1, wherein, The identification of the wet and dry states of the mesh in step S1 comprises the following steps: introducing a critical water depth threshold as a judgment variable, when the water depth in the mesh is greater than or equal to the critical threshold, the mesh is determined as a wet mesh, otherwise, the mesh is determined as a dry mesh; the interface flux processing comprises the following cases: (1) both sides of the mesh are dry meshes, or the liquid level of the wet mesh is lower than that of the dry mesh, and there is no flux at the interface; (2) both sides of the mesh are wet meshes, or the liquid level of the wet mesh is greater than or equal to that of the dry mesh, and there is flux at the interface.

3. The method of claim 1, wherein, The adaptation of the complex boundary conditions in step S1 comprises the following steps: when all the three edges of the mesh are no-flux boundaries, the mesh does not participate in the calculation; when all the three edges of the mesh are flux boundaries, the mesh is calculated normally; when part of the three edges of the mesh are flux boundaries, the mesh only calculates the mass equation and does not calculate the momentum equation.

4. The method of claim 1, wherein, The discretization of the finite volume method in step S2 comprises the following steps: the volume integral of the control equation is converted into a boundary integral by using the Gauss divergence theorem, and the flow update formula and the water depth update formula where Q i is the flow vector of cell i; A i is the area of cell i; F conv , F pres , F diff are the fluxes of convection, pressure and diffusion terms, respectively; Δt is the time step; h i is the water depth of cell i, Q f is the flow of the boundary edge, n f is the edge normal vector, L f is the edge length.

5. The method of claim 1, wherein, are obtained after discretization.

6. The method of claim 1, wherein, The dynamic adjustment of the time step in step S4 is based on the CFL condition, combined with the local wave speed and the minimum characteristic length of the mesh element to dynamically calculate the maximum allowed time step, the time step is automatically reduced when the flow velocity gradient is large or the mesh element is small, and vice versa. In step S2, the NumPy vector operation and Numba library just-in-time compilation technology are used to optimize the core calculation process, the cell edge value calculation and flux accumulation nested loop logic are rewritten into matrix operations, and the @numba.jit decorator is added to the gradient reconstruction and time step iteration core functions for machine code compilation.

7. The method of claim 1, wherein, The step S1 further comprises an adaptive grid refinement function, which identifies high variation areas in the flow field based on the length of the flow velocity gradient or the length of the water depth gradient, and subdivides the grid cells whose gradient values exceed a set threshold, so as to divide the triangular cells into four sub-triangular cells.

8. An explicit unstructured two-dimensional hydrodynamic model construction system, characterized by, The system is based on the explicit unstructured two-dimensional hydrodynamic model construction method of claim 1, comprising: a grid and data processing module for constructing an unstructured grid, identifying the wet and dry states of the grid based on a critical water depth threshold, and adapting to complex boundaries according to the flux characteristics of the boundary edges; a model construction module for realizing flow field prediction, correction and water depth update based on the conservation type two-dimensional shallow water equation, the finite volume method for discretizing the control equation, and the explicit pressure-velocity coupling mechanism; a parameter configuration module for setting boundary conditions, initial conditions and dynamic time step parameters; a calculation iteration module for dynamic iterative calculation using an explicit time stepping method; a result output module for outputting simulation results and performing visual processing and flow motion characteristic analysis.

Images