Surface large-scale hydrodynamic simulation method based on shallow water equation set

Abstract

The invention discloses an earth surface large-scale hydrodynamic simulation method based on a shallow water equation set, and belongs to the technical field of hydromechanics, and the method comprises the steps: introducing a bed bottom source item, and constructing a shallow water lattice Boltzmann equation; constructing a dry-wet boundary treatment equation in combination with the single relaxation model; constructing a model to be simulated, and setting a computational domain according to the constructed model; setting simulation parameters of the to-be-simulated model; setting boundary conditions of the to-be-simulated model; and simulating the dynamic phenomenon of the to-be-simulated model according to the constructed shallow water lattice Boltzmann equation and dry-wet boundary treatment equation. By the adoption of the method, on the premise that mass-momentum conservation is guaranteed, the calculation stability of dynamic boundary scenes such as the dam break wave front edge and intertidal zone submerging is remarkably improved, and calculation time is shortened.

Description

Technical Field

[0001] The present invention relates to the technical field of fluid mechanics, and in particular to a surface large-scale hydrodynamic simulation method based on a shallow water equation group. Background Art

[0002] The shallow water equations (SWEs), as the core mathematical model for describing surface hydrodynamic processes, originated from the fundamental research on fluid mechanics by Euler and Lagrange in the 18th century. This set of equations is obtained by vertically integrating the Navier-Stokes equations and is applicable to flows where the horizontal scale is much larger than the vertical scale. The classical shallow water equations consist of the mass conservation equation and the momentum conservation equation, which can effectively characterize the mass transfer, momentum exchange, and energy conversion processes in large-scale surface water flows. Since the mid-20th century, with the rise of computational fluid dynamics, the shallow water equations have demonstrated unique advantages in large-scale hydrodynamic simulations such as flood evolution, estuary dynamics, and tsunami propagation, becoming an important tool in the fields of hydrology, oceanography, and environmental engineering.

[0003] Large-scale surface hydrodynamic simulations face unique challenges such as spatial heterogeneity and multi-scale coupling. Traditional numerical solutions to shallow water equations (such as the finite volume method and the finite difference method) face significant technical bottlenecks. First, in the treatment of bed source terms in areas with sudden changes in topography, traditional methods often use simplified discretization schemes, resulting in inaccurate nonlinear coupling between the topographic gradient term and the momentum equation. For example, in complex terrain areas such as steep slopes and stepped riverbeds, existing algorithms are unable to accurately represent the characteristics of sudden changes in elevation, which can easily cause non-physical oscillations in the velocity field, resulting in decreased computational stability and loss of accuracy. Second, traditional dry-wet boundary treatment methods use fixed criteria, which makes it difficult to dynamically capture the migration of the water-land interface. Especially in the rapid evolution stage of floods, non-physical diffusion phenomena are prone to occur at the submerged boundary, significantly affecting the reliability of the prediction of the instantaneous submerged range. In addition, in multi-physics field coupling calculations, traditional methods need to rely on iterative strategies, which not only greatly increases the computational complexity, but also deteriorates the mass conservation due to improper matching of the relaxation mechanism. Summary of the Invention

[0004] The purpose of the present invention is to provide a surface large-scale hydrodynamic simulation method based on a shallow water equation group to solve the technical problems in the background technology.

[0005] To achieve the above object, the present invention provides a method for simulating large-scale surface hydrodynamics based on a shallow water equation system, the steps comprising:

[0006] S1. Introduce the bed source term and construct the shallow water lattice Boltzmann equation;

[0007] S2. Construct dry-wet boundary treatment equations in combination with the single relaxation model;

[0008] S3. Construct a model to be simulated and set a computational domain according to the constructed model;

[0009] S4, setting simulation parameters of the model to be simulated;

[0010] S5. Setting the boundary conditions of the model to be simulated;

[0011] S6. The dynamic phenomena of the model to be simulated are simulated based on the constructed shallow water lattice Boltzmann equation and the dry-wet boundary treatment equation.

[0012] Preferably, step S1 specifically includes:

[0013] Construct the shallow water lattice Boltzmann equation, the formula is:

[0014]

[0015] Where, f α is the particle distribution function; x is a space vector, in two-dimensional space x = (x, y); t is time; Δt is the time step; e α is the particle velocity vector; e αj It is e α The component in the jth direction; e = Δx / Δt, Δx is the grid size; z b is the bed elevation; τ is the single relaxation time; is the local equilibrium distribution function; ω α represents the weighting coefficient;

[0016]

[0017] Where, F j is the component of the volume force acting on the fluid in the j direction, h is the water depth, which is the height from the bed to the water surface; u i is the depth-averaged velocity; τ ωj is the wind shear stress in the jth direction; τ bj is the bed shear stress; ρ is the water density; u j is the component of the depth-averaged velocity in the jth direction; u y is the component of the depth-averaged velocity in the y direction; u x is the component of the depth-averaged velocity in the x direction; C b is the bed friction coefficient; Ω is the Coriolis parameter; δ ij is the Kronecker function, and its formula is:

[0018]

[0019] Redefine the general local equilibrium distribution function as follows:

[0020]

[0021] Where A and B are unknown constants; g is the acceleration due to gravity; e αi is the particle velocity vector e α The component in the i direction (usually the x or y direction).

[0022] Preferably, step S2 specifically includes:

[0023] Combined with the single relaxation model to construct the dry-wet boundary treatment equation, in the shallow water lattice Boltzmann equation, the flow passes through the distribution function f α The evolution of the collision and migration process of (x, t) is as follows:

[0024]

[0025] Where, F α is the external force term;

[0026] Decompose the distribution function into equilibrium and non-equilibrium disturbance terms, the formula is:

[0027]

[0028] Where, ∈ is a dimensionless small parameter; at the dry unit h=0, velocity u=0, is an unbalanced term, which is transmitted through spatial gradient;

[0029] Discretize time and space, perform Taylor expansion on time and space, and approximate partial derivatives in discrete format. The formula is:

[0030]

[0031] Where, Represents the distribution function of the wet unit. At the wet unit, h>0; represents the distribution function of the stem unit;

[0032] Calculate the wet cell distribution function The spatial gradient of Spatial ladder calculation of dry unit distribution function

[0033] Preferably, the calculation of the wet unit distribution function The spatial gradient of the stem cell distribution function is reversed The spatial gradients include:

[0034] Calculate the wet cell distribution function The formula for the spatial gradient is:

[0035]

[0036] Where α is the discrete velocity direction index; Δx is the spatial step size, is the distribution function value on the right side of the right wet unit, is the distribution function value on the left side of the wet unit;

[0037] Dry cell distribution function The formula after spatial gradient correction is:

[0038]

[0039] Where, is the correction factor, is the gradient information.

[0040] Preferably, the model to be simulated in step S3 is a dam break wave scouring cone group and a slope revetment model, including upstream and downstream areas and a downstream slope area, the downstream area includes the cone group, and the calculation domain includes two variables, length and width.

[0041] Preferably, the simulation parameters in step S4 include calculation time step, fluid density, relaxation factor, gravitational acceleration, Manning coefficient and wind speed.

[0042] Preferably, the upstream region and the left and right sides of the flow field are set as no-slip solid wall boundaries, and the cone group and the downstream slope region are set as non-equilibrium extrapolation boundary conditions.

[0043] Therefore, the present invention adopts the above-mentioned large-scale surface hydrodynamic simulation method based on the shallow water equations, which has the following beneficial effects:

[0044] (1) An improved bed source term is introduced into the original shallow water lattice Boltzmann equation model framework, and a dry-wet boundary treatment scheme is constructed in combination with a single relaxation model. By automatically determining the dry-wet boundary conditions, the distribution function evolution can adaptively track the migration of the water-land interface;

[0045] (2) The distribution function of the dry-wet transition zone is dynamically reconstructed based on gradient information, and the spatial gradient of the adjacent wet unit is used to correct the state of the dry unit. This not only avoids the non-physical errors introduced by artificial parameters, but also effectively suppresses the numerical singularity caused by the water depth approaching zero. At the same time, an adaptive discretization format based on the particle velocity vector is introduced, and the numerical representation ability of the terrain mutation is enhanced through the directional projection difference operator, which significantly reduces the momentum flux calculation error in the steep slope and stepped terrain area.

[0046] (3) A single relaxation model is used to decouple the time scale differences of physical processes and achieve efficient synchronous evolution of exogenous terms through feature space separation, maintaining mass conservation while ensuring the advantage of a large time step;

[0047] (4) The local dependence characteristics of the distribution function are highly compatible with the natural parallel architecture of the shallow water lattice Boltzmann method. While ensuring the conservation of mass-momentum, it significantly improves the computational stability of dynamic boundary scenarios such as the dam break wave front and intertidal zone inundation.

[0048] The technical solution of the present invention is further described in detail below through the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 This is a flow chart of the method of Example 1 of the present invention;

[0050] Figure 2 This is a diagram of a dam-break wave scouring cone group and a slope bank protection model according to Example 1 of the present invention;

[0051] Figure 3 This is a simulation process diagram of the Nth calculation time step in Example 1 of the present invention;

[0052] Figure 4 The water depth and flow field change diagram at four moments when the upstream and downstream water depth is 0.3m in Example 1 of the present invention;

[0053] Figure 5 Graphs showing changes in water depth and flow field at four moments when the downstream water depth is 0.1 m according to Example 1 of the present invention;

[0054] Figure 6 Schematic diagram of the flow field at four different moments when wind blows on a circular lake surface, as simulated in Example 2 of the present invention;

[0055] Figure 7 This is a schematic diagram of the bottom slope height of a simulated river section according to Example 3 of the present invention;

[0056] Figure 8 Schematic diagram of the flow field at four different moments of the dam-break flood dynamic simulation model in Example 3 of the present invention. DETAILED DESCRIPTION

[0057] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings 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, and therefore cannot be understood as limiting the present invention.

[0058] Example 1

[0059] Reference Figure 1 The present invention provides a method for simulating large-scale surface hydrodynamics based on a shallow water equation group. The method is described by taking the hydrodynamic calculation of multi-cone scour of a dam-break wave as an example. The steps include:

[0060] S1. Introduce the bed source term and construct the shallow water lattice Boltzmann equation.

[0061] Specifically, the shallow water lattice Boltzmann equation is:

[0062]

[0063] Where, f α is the particle distribution function; x is a space vector, in two-dimensional space x = (x, y); t is time; Δt is the time step; e α is the particle velocity vector; e αj It is e α The component in the jth direction; e = Δx / Δt, Δx is the grid size; z b is the bed elevation; τ is the single relaxation time; is the local equilibrium distribution function; ω α represents the weighting coefficient;

[0064]

[0065] Where, F j is the component of the volume force acting on the fluid in the j direction, h is the water depth, which is the height from the bed to the water surface; u i is the depth-averaged velocity; τ ωj is the wind shear stress in the jth direction; τ bj is the bed shear stress; ρ is the water density u j is the component of the depth-averaged velocity in the jth direction; u y is the component of the depth-averaged velocity in the y direction; u x is the component of the depth-averaged velocity in the x direction; C b is the bed friction coefficient; Ω is the Coriolis parameter; δ ij is the Kronecker function, and its formula is:

[0066]

[0067] Redefine the general local equilibrium distribution function as follows:

[0068]

[0069] Where A and B are unknown constants; g is the acceleration due to gravity, g = 9.81 m / s 2 ;e αi is the particle velocity vector e α The component in the i direction (usually the x or y direction); A and B can also be expressed as λ α Expressed as follows:

[0070]

[0071] The physical variables defining water depth and velocity are:

[0072] h=∑ a f a ;

[0073] u i h=∑ a e ai f a .

[0074] S2. Combine the single relaxation model to construct the dry-wet boundary treatment equation.

[0075] Specifically, in the shallow water lattice Boltzmann equation, the flow is distributed through the distribution function f α The evolution of the collision and migration process of (x, t) is as follows:

[0076]

[0077] Where, F α is an external force term, which can be bottom friction or terrain source term, etc.

[0078] The distribution function of the dry unit (water depth h = 0) cannot be directly obtained through the equilibrium state The calculation requires the reconstruction of the gradient information of the adjacent wet cells. The distribution function is decomposed into equilibrium and non-equilibrium disturbance terms, and the formula is:

[0079]

[0080] Where, ∈ is a dimensionless small parameter; at the dry unit h=0, velocity u=0, It is an unbalanced term and needs to be transferred through spatial gradient;

[0081] Discretize time and space, perform Taylor expansion on time and space, and approximate partial derivatives in discrete format. The formula is:

[0082]

[0083] Where, Represents the distribution function of the wet unit. At the wet unit, h>0; represents the distribution function of the stem unit;

[0084] Calculate the wet cell distribution function The spatial gradient of Spatial ladder calculation of dry unit distribution function

[0085] In wet cells (h>0), calculate the wet cell distribution function The formula for the spatial gradient is:

[0086]

[0087] Where α is the discrete velocity direction index; Δx is the spatial step size, is the distribution function value on the right side of the right wet unit, is the distribution function value on the left side of the wet unit;

[0088] Dry cell distribution function The formula after spatial gradient correction is:

[0089]

[0090] In the formula, the first term on the right side of the equation is the average value (weighted average) of the distribution function of adjacent wet units, and the second term is the gradient correction term. is the correction factor, is the gradient information.

[0091] The wet-dry boundary treatment method significantly optimizes the boundary condition adaptability of the shallow water lattice Boltzmann method through a physically driven mechanism. The wet-dry boundary treatment dynamically reconstructs the distribution function of the dry-wet transition zone based on gradient information and modifies the state of dry cells using the spatial gradients of adjacent wet cells. This avoids non-physical errors introduced by artificial parameters and effectively suppresses numerical singularities caused by water depth tending to zero. Its local dependence is highly compatible with the inherently parallel architecture of the shallow water lattice Boltzmann method. While ensuring mass-momentum conservation, it significantly improves the computational stability and simulation accuracy for dynamic boundary scenarios such as dam-break wave fronts and intertidal zone inundation.

[0092] S3. Construct the model to be simulated and set the calculation domain according to the constructed model.

[0093] Specifically, a dam-break wave scouring cone group and slope revetment model are constructed, such as Figure 2 The aforementioned method includes upstream and downstream areas and a downstream slope area, wherein the initial water depth upstream is 0.3 m, the initial water depth downstream is 0.1 m, the downstream area includes a cone group, and the calculation domain is set according to the constructed model. The length L in the calculation domain is 600 and the width W is 200, both of which are grid units.

[0094] S4. Set simulation parameters of the model to be simulated.

[0095] Specifically, the specific simulation parameters of the dam break wave scour cone group and slope revetment model are set. The simulation calculation time step, fluid density, relaxation factor, gravitational acceleration, Manning coefficient and wind speed can be set by variables such as t, rho, tau, gravity, manning, and windspeed. Among them, t is 20,000 time steps, rho is 1, tau is 0.53, gravity is 9.81, manning is 0.04, and windspeed is 0. It should be noted that if high-precision large-scale hydrodynamic phenomena are to be simulated, the settings of these variables need to be converted according to certain unit conversion principles based on real physical quantities.

[0096] S5. Set the boundary conditions of the model to be simulated.

[0097] Specifically, the upstream area and the left and right sides of the flow field of the dam break wave scouring cone group and slope revetment model are set as no-slip solid wall boundaries, and the cone group and downstream slope area are set as non-equilibrium extrapolation boundary conditions.

[0098] S6. Based on the constructed shallow water lattice Boltzmann equation and the dry-wet boundary treatment equation, the dynamic phenomena of the dam-break wave scouring the cone group and the slope revetment model are simulated. The Nth calculation time step of the simulation is as follows: Figure 3 As shown, the calculation is stopped when the flow field velocity change is less than the judgment condition, that is, when ΔU ≤ 0.01, the simulation results are output.

[0099] Specifically, the water depth and flow field changes at four moments, t = 2000, t = 6000, t = 10000, and t = 16000, were simulated when the upstream water depth was 0.3m and the upstream water depth was 0.3m and the downstream water depth was 0.1m, as shown in the figure below. Figure 4 and Figure 5 As shown in the figure, at t = 2000, the water flows from upstream to downstream. At this time, the wave front does not completely contact the large cone. As the flow continues, the dam-break wave bypasses the cone and, at t = 6000, bypasses the two smaller cones, further changing the flow field. At t = 10000, the dam-break wave first evolves to the rightmost slope, forming a reflected wave. This phenomenon can be observed based on the streamlines. This shows that the method of the present invention improves the computational stability of dynamic boundary scenarios such as the dam-break wave front and intertidal zone inundation.

[0100] Example 2

[0101] The constructed shallow water lattice Boltzmann equation and the dry-wet boundary treatment equation are used to simulate the wind blowing on the circular lake surface. A computational domain of 200 grids long and 200 grids wide is set. The wind speed in the x-direction is set to 0 and the wind speed in the y-direction is set to 9 m / s. Assuming that it is not affected by the rotation of the earth, the flow field diagrams at different times are shown as follows: Figure 6 As shown, through Figure 6, it can be seen that the simulation equation constructed by the present invention can accurately capture the changes in water surface flow velocity. At T = 500, the water surface streamlines are upward as the wind blows; at T = 1500, the streamlines turn downward, forming clockwise and counterclockwise vortices on the left and right sides of the lake, respectively. As the wind continues to blow, the vortices grow larger and alternately change the direction of the streamlines; at T = 15000, two large vortices form on each side, and the streamline direction no longer changes with the wind. These simulation results are consistent with real-world physical laws, demonstrating the effectiveness of the present method in simulating large-scale surface hydrodynamic models.

[0102] Example 3

[0103] A dam-break flood simulation was conducted on a section of the Minjiang River. A dynamic simulation model of the dam-break flood was constructed, and the flow rate and inundation area of ​​the dam-break flood at different times were analyzed. The bottom slope height diagram of the section of the dam-break flood dynamic simulation model is shown in the figure below. Figure 7 shown.

[0104] By solving the shallow water lattice Boltzmann equations, importing the measured data of the river section, and using adaptive grid technology to achieve a refined expression of the terrain features. In the hypothetical dam breach scenario of a cascade power station downstream of the Zipingpu Reservoir, the breach width is set to 200 meters and the initial flow peak is 9000m 3 / s, and conduct dynamic simulation of flood evolution process for up to 1 hour. Figure 8 As shown, a flood velocity map at different times is displayed, including the flooded area.

[0105] Simulation results show that the dam-break wave generated impact velocities of up to 5 m / s in narrow river sections, flooding some low-lying areas within 40 minutes of the dam breach, with water reaching a maximum depth of 4.5 meters. By outputting parameters such as flooding depth, velocity vector field, and flood arrival time, the interaction between the dam-break flood and topography was analyzed.

[0106] The method in this paper takes 23 minutes to run on a computer equipped with an Intel i7-12700 processor, while the well-known flood simulation software Mike 21 takes 51 minutes on the same platform. Engineering applications have shown that this method significantly reduces computational time in flood simulations in complex terrain and improves the reliability of predictions of key parameters such as flood depth and flow velocity distribution, providing more accurate technical support for disaster prevention and control decisions.

[0107] Therefore, the present invention adopts the above-mentioned large-scale surface hydrodynamic simulation method based on the shallow water equations, which significantly improves the calculation stability of dynamic boundary scenes such as the dam break wave front and intertidal zone flooding, and shortens the calculation time while ensuring the conservation of mass-momentum.

[0108] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention rather than to limit the same. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that they can still modify or replace the technical solutions of the present invention with equivalents, and these modifications or equivalent replacements cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A method for simulating large-scale surface hydrodynamics based on shallow water equations, characterized in that the steps include: S1. Introduce the bed source term and construct the shallow water lattice Boltzmann equation; S2. Construct dry-wet boundary treatment equations in combination with the single relaxation model; S3. Construct a model to be simulated and set a computational domain according to the constructed model; S4, setting simulation parameters of the model to be simulated; S5. Setting the boundary conditions of the model to be simulated; S6. The dynamic phenomena of the model to be simulated are simulated based on the constructed shallow water lattice Boltzmann equation and the dry-wet boundary treatment equation.

2. The surface large-scale hydrodynamic simulation method based on shallow water equations according to claim 1 is characterized in that: Step S1 specifically includes: Construct the shallow water lattice Boltzmann equation, the formula is: Where, f α is the particle distribution function; x is a space vector, in two-dimensional space x = (x, y); t is time; Δt is the time step; e α is the particle velocity vector; e αj It is e α The component in the jth direction; e = Δx / Δt, Δx is the grid size; z b is the bed elevation; τ is the single relaxation time; is the local equilibrium distribution function; ω α represents the weighting coefficient; Where, F j is the component of the volume force acting on the fluid in the j direction, h is the water depth, which is the height from the bed to the water surface; u i is the depth-averaged velocity; τ ωj is the wind shear stress in the jth direction; τ bj is the bed shear stress; ρ is the water density; u j is the component of the depth-averaged velocity in the jth direction; u y is the component of the depth-averaged velocity in the y direction; u x is the component of the depth-averaged velocity in the x direction; C b is the bed friction coefficient; Ω is the Coriolis parameter; δ ij is the Kronecker function, and its formula is: Redefine the general local equilibrium distribution function as follows: Where A and B are unknown constants; g is the acceleration due to gravity; e αi is the particle velocity vector e α Component in the i direction.

3. The surface large-scale hydrodynamic simulation method based on shallow water equations according to claim 2 is characterized in that: Step S2 specifically includes: Combined with the single relaxation model to construct the dry-wet boundary treatment equation, in the shallow water lattice Boltzmann equation, the flow passes through the distribution function f α The evolution of the collision and migration process of (x, t) is as follows: Where, F α is the external force term; Decompose the distribution function into equilibrium and non-equilibrium disturbance terms, the formula is: Where, ∈ is a dimensionless small parameter; at the dry unit h=0, velocity u=0, is an unbalanced term, which is transmitted through spatial gradient; Discretize time and space, perform Taylor expansion on time and space, and approximate partial derivatives in discrete format. The formula is: Where, Represents the distribution function of the wet unit. At the wet unit, h>0; represents the distribution function of the stem unit; Calculate the wet cell distribution function The spatial gradient of Spatial ladder calculation of dry unit distribution function 4. The surface large-scale hydrodynamic simulation method based on shallow water equations according to claim 3 is characterized in that: The calculated wet cell distribution function The spatial gradient of the stem cell distribution function is reversed The spatial gradients include: Calculate the wet cell distribution function The formula for the spatial gradient is: Where α is the discrete velocity direction index; Δx is the spatial step size, is the distribution function value on the right side of the right wet unit, is the distribution function value on the left side of the wet unit; Dry cell distribution function The formula after spatial gradient correction is: Where, is the correction factor, is the gradient information.

5. The method for simulating large-scale surface hydrodynamics based on shallow water equations according to claim 1, characterized in that: The model to be simulated in step S3 is a dam break wave scouring cone group and a slope revetment model, including upstream and downstream areas and a downstream slope area. The downstream area includes the cone group, and the calculation domain includes two variables: length and width.

6. The method for simulating large-scale surface hydrodynamics based on shallow water equations according to claim 1, characterized in that: The simulation parameters in step S4 include calculation time step, fluid density, relaxation factor, gravitational acceleration, Manning coefficient and wind speed.

7. The surface large-scale hydrodynamic simulation method based on shallow water equations according to claim 1 is characterized in that: The boundary conditions in step S5 include: the upstream area and the left and right sides of the flow field are set as no-slip solid wall boundaries, and the cone group and the downstream slope area are set as non-equilibrium extrapolation boundary conditions.

Images