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

Abstract

The application discloses a kind of ground surface large-scale hydrodynamics simulation simulation methods based on shallow water equation group, belong to the technical field of fluid mechanics, steps include: introducing bed bottom source term, shallow water lattice Boltzmann equation is constructed;Combining single relaxation model constructs dry-wet boundary processing equation;The model to be simulated is constructed, and the calculation domain is set according to the constructed model;Simulation parameters of the model to be simulated are set;The boundary condition of the model to be simulated is set;According to the shallow water lattice Boltzmann equation and dry-wet boundary processing equation constructed, the dynamic phenomenon of the model to be simulated is simulated.The application uses the above method, on the premise of guaranteeing mass-momentum conservation, significantly improves the calculation stability of dam-break wave front, intertidal zone submergence and other dynamic boundary scenarios, and shortens the time consumption of calculation.

Description

Technical Field

[0001] This invention relates to the field of fluid mechanics, and in particular to a method for large-scale hydrodynamic simulation of the Earth's surface based on shallow water equations. Background Technology

[0002] Shallow water equations (SWEs), as core mathematical models describing surface hydrodynamic processes, originated from the fundamental research on fluid mechanics by Euler and Lagrange in the 18th century. These equations are obtained through vertical integration of the Navier-Stokes equations and are applicable to flows with horizontal scales much larger than vertical scales. The classical shallow water equations consist of mass conservation and momentum conservation equations, effectively characterizing mass transport, momentum exchange, and energy conversion processes in large-scale surface water flows. Since the mid-20th century, with the rise of computational fluid dynamics, shallow water equations have demonstrated unique advantages in large-scale hydrodynamic simulations such as flood evolution, estuarine dynamics, and tsunami propagation, becoming an important tool in hydrology, oceanography, and environmental engineering.

[0003] Large-scale hydrodynamic simulations of the Earth's surface face unique challenges such as spatial heterogeneity and multi-physics-scale coupling, and traditional numerical methods for solving shallow water equations (such as the finite volume method and the finite difference method) encounter significant technical bottlenecks. First, in handling the source term at the bottom of the bed in areas with abrupt topographic changes, traditional methods often employ simplified discretization schemes, leading to inaccurate nonlinear coupling between the topographic gradient term and the momentum equation. For example, in complex topographic regions such as steep slopes and stepped riverbeds, existing algorithms, unable to accurately characterize elevation changes, are prone to non-physical oscillations in the velocity field, resulting in decreased computational stability and accuracy loss. Second, traditional methods for handling wet and dry boundaries use fixed criteria, making it difficult to dynamically capture the migration of the land-water interface. Especially during the rapid evolution phase of floods, non-physical diffusion phenomena easily occur at the inundation boundary, significantly affecting the reliability of instantaneous inundation range predictions. Furthermore, in multi-physics coupled calculations, traditional methods rely on iterative strategies, which not only significantly increases computational complexity but also leads to deterioration of mass conservation due to improper matching of relaxation mechanisms. Summary of the Invention

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

[0005] To achieve the above objectives, this invention provides a method for large-scale surface hydrodynamic simulation based on shallow water equations, comprising the following steps:

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

[0007] S2. Construct equations for handling wet and dry boundaries by combining a single relaxation model;

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

[0009] S4. Set the simulation parameters for the model to be simulated;

[0010] S5. Set the boundary conditions for the model to be simulated;

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

[0012] Preferably, step S1 specifically includes:

[0013] The Boltzmann equations for shallow-water lattice systems are constructed as follows:

[0014]

[0015] In the formula, f α It is the particle distribution function; x is a spatial vector, in two-dimensional space x = (x, y); t is time; Δt is the time step; e α It is the particle velocity vector; e αj It is e α The component in the j-th direction; e = Δx / Δt, where Δx is the lattice size; z b τ is the bed surface elevation; τ is the single relaxation time. It is a local equilibrium distribution function; ω α Indicates the weighting coefficient;

[0016]

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

[0018]

[0019] The general local equilibrium distribution function is redefined as follows:

[0020]

[0021] In the formula, A and B are unknown constants; g is the acceleration due to gravity; e αi It 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] Combining a single relaxation model to construct equations for treating wet and dry boundaries, in the shallow water lattice Boltzmann equations, the flow passes through the distribution function f α The evolution of the collision and migration process at (x,t) is expressed by the following formula:

[0024]

[0025] In the formula, F α External force term;

[0026] The distribution function is decomposed into equilibrium and non-equilibrium disturbance terms, as shown in the formula:

[0027]

[0028] In the formula, ∈ represents a dimensionless small parameter; at the dry element, h = 0 and velocity u = 0. This is a non-equilibrium term, propagated through spatial gradient;

[0029] Discretize the spacetime, perform Taylor expansion on time and space, and combine the discretization scheme to approximate the partial derivatives. The formula is as follows:

[0030]

[0031] In the formula, Let h represent the distribution function of the wet element, where h > 0 at the wet element. Represents the distribution function of the dry unit;

[0032] Calculate the distribution function of the wetted element The spatial gradient, based on the wet unit distribution function Spatial ladder calculation dry element distribution function

[0033] Preferably, the calculation of the wet unit distribution function Spatial gradient, inversely deriving the dry unit distribution function Spatial gradients include:

[0034] Calculate the distribution function of the wetted element The formula for the spatial gradient is:

[0035]

[0036] In the formula, α is the index of the discrete velocity direction; Δx is the spatial step size. The distribution function value is the value on the right side of the right wetted element. The value of the distribution function on the left side of the wet cell;

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

[0038]

[0039] In the formula, For correction factor, This represents gradient information.

[0040] Preferably, the model to be simulated in step S3 is a dam break wave scour cone group and a slope revetment model, including upstream, downstream and downstream slope regions, the downstream region includes the cone group, and the computational domain includes two variables: length and width.

[0041] Preferably, the simulation parameters in step S4 include the 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 non-slip solid wall boundaries, and the conical group and the downstream slope region are set as non-equilibrium extrapolation boundary conditions.

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

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

[0045] (2) Based on gradient information, the distribution function of the dry-wet transition zone is dynamically reconstructed. The state of the dry unit is corrected by the spatial gradient of the neighboring wet unit. This avoids the non-physical error introduced by artificial parameters and effectively suppresses the numerical singularity caused by the water depth approaching zero. At the same time, an adaptive discrete format based on particle velocity vector is introduced. The numerical representation ability of terrain change is enhanced by the directional projection difference operator, which significantly reduces the momentum flux calculation error in steep slope and stepped terrain areas.

[0046] (3) A single relaxation model is used to decouple the time scale differences of the physical process. The efficient synchronous evolution of the exogenous term is achieved through feature space separation, which maintains the mass conservation while ensuring the advantage of large time step.

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

[0048] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0049] Figure 1 This is a flowchart of the method in Embodiment 1 of the present invention;

[0050] Figure 2 This is a model diagram of the dam break wave scour cone group and slope revetment in Embodiment 1 of the present invention;

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

[0052] Figure 4 The diagram shows the changes in water depth and flow field at four times when the upstream and downstream water depth is 0.3m in Embodiment 1 of the present invention.

[0053] Figure 5 This is a diagram showing the changes in water depth and flow field at four times when the downstream water depth is 0.1m in Embodiment 1 of the present invention.

[0054] Figure 6 This is a schematic diagram of the flow field at four different moments when wind blows across a circular lake surface, as shown in Embodiment 2 of the present invention.

[0055] Figure 7 This is a schematic diagram of the simulated riverbed slope height in Embodiment 3 of the present invention;

[0056] Figure 8 This is a schematic diagram of the flow field at four different times in the dynamic simulation model of dam-break flood in Embodiment 3 of the present invention. Detailed Implementation

[0057] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments, and therefore should not be construed as limiting the present invention.

[0058] Example 1

[0059] Reference Figure 1 This invention provides a large-scale surface hydrodynamic simulation method based on shallow water equations, using the hydrodynamic calculation of dam-break wave multi-cone scour as an example. The steps include:

[0060] S1. Introduce the bottom source term to construct the shallow water lattice Boltzmann equation.

[0061] Specifically, the formula for the Boltzmann equation in shallow water lattice is as follows:

[0062]

[0063] In the formula, f α It is the particle distribution function; x is a spatial vector, in two-dimensional space x = (x, y); t is time; Δt is the time step; e α It is the particle velocity vector; e αj It is e α The component in the j-th direction; e = Δx / Δt, where Δx is the lattice size; z b τ is the bed surface elevation; τ is the single relaxation time. It is a local equilibrium distribution function; ω α Indicates the weighting coefficient;

[0064]

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

[0066]

[0067] The general local equilibrium distribution function is redefined as follows:

[0068]

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

[0070]

[0071] The physical variables of water depth and velocity are defined by the following formulas:

[0072] h = ∑ a f a ;

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

[0074] S2. Construct equations for handling wet and dry boundaries by combining a single relaxation model.

[0075] Specifically, in the shallow water lattice Boltzmann equation, the flow passes through the distribution function f α The evolution of the collision and migration process at (x,t) is expressed by the following formula:

[0076]

[0077] In the formula, F α This is an external force term, which can be bottom friction or topographic source term, etc.

[0078] The distribution function of the dry element (water depth h = 0) cannot be directly obtained through equilibrium state. The calculation requires reconstruction using gradient information from neighboring wetted cells. The distribution function is decomposed into equilibrium and non-equilibrium perturbation terms, as shown in the formula:

[0079]

[0080] In the formula, ∈ represents a dimensionless small parameter; at the dry element, h = 0 and velocity u = 0. Since it is a non-equilibrium term, it needs to be propagated through spatial gradient;

[0081] Discretize the spacetime, perform Taylor expansion on time and space, and combine the discretization scheme to approximate the partial derivatives. The formula is as follows:

[0082]

[0083] In the formula, Let h represent the distribution function of the wet element, where h > 0 at the wet element. Represents the distribution function of the dry unit;

[0084] Calculate the distribution function of the wetted element The spatial gradient, based on the wet unit distribution function Spatial ladder calculation dry element distribution function

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

[0086]

[0087] In the formula, α is the index of the discrete velocity direction; Δx is the spatial step size. The distribution function value is the value on the right side of the right wetted element. The value of the distribution function on the left side of the wet cell;

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

[0089]

[0090] In the formula, the first term on the right-hand side is the average (weighted average) of the distribution functions of adjacent wetted units, and the second term is the gradient correction term. For correction factor, This represents gradient information.

[0091] The wet-dry boundary treatment method significantly optimizes the boundary condition adaptability of the shallow-water lattice Boltzmann method through a physical driving mechanism. Based on gradient information, the wet-dry boundary treatment dynamically reconstructs the distribution function of the wet-dry transition zone and uses the spatial gradient of neighboring wet cells to correct the state of dry cells. This avoids non-physical errors introduced by artificial parameters and effectively suppresses numerical singularities caused by water depth approaching zero. Its local dependency characteristics are highly compatible with the natural parallel architecture of the shallow-water lattice Boltzmann method. While ensuring mass-momentum conservation, it significantly improves the computational stability of dynamic boundary scenarios such as dam-break wavefronts and intertidal inundation, and also enhances simulation accuracy.

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

[0093] Specifically, models of dam-break wave scour cone groups and sloping revetments are constructed, such as... Figure 2 The system includes an upstream region, a downstream region, and a downstream slope region. The initial water depth in the upstream region is 0.3m, and the initial water depth in the downstream region is 0.1m. The downstream region includes a group of cones. The computational domain is set according to the constructed model. The length of the computational domain is L=600 and the width is W=200, both of which are grid units.

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

[0095] Specifically, the simulation parameters for the dam-break wave scouring 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. Here, t is 20000 time steps, rho is 1, tau is 0.53, gravity is 9.81, manning is 0.04, and windspeed is 0. It is important to note that to simulate high-precision, large-scale hydrodynamic phenomena, these variables need to be converted according to certain unit conversion principles based on the actual physical quantities.

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

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

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

[0099] Specifically, simulations were performed at t=2000, t=6000, t=10000, and t=16000, with upstream water depth of 0.3m and downstream water depth of 0.1m, respectively, showing the changes in water depth and flow field at four time points: upstream water depth of 0.3m and downstream water depth of 0.1m. (See diagrams for reference.) Figure 4 and Figure 5 As shown, at t=2000, the water flows from upstream to downstream. At this time, the wave front has not completely contacted the large cone. As the flow continues, the dam-break wave bypasses the cone, and at t=6000, it bypasses two smaller cones, further changing the flow field. At t=10000, the dam-break wave evolves to the rightmost slope for the first time, forming a reflected wave. This phenomenon can be observed based on the streamlines. Therefore, it can be concluded that the method of this invention improves the computational stability of dynamic boundary scenarios such as the dam-break wave front and intertidal submersion.

[0100] Example 2

[0101] The constructed shallow-water lattice Boltzmann equation and wet / dry boundary treatment equations were used to simulate wind blowing across a circular lake surface. A computational domain of 200 grids in length and 200 grids in width was defined, with the wind speed in the x-direction set to 0 and the wind speed in the y-direction set to 9 m / s. Assuming no influence from the Earth's rotation, the flow field diagrams at different times are shown below. Figure 6 As shown, through Figure 6As can be seen, the simulation equations constructed in this invention can accurately capture changes in water surface velocity. At T=500, the streamlines move upward with the wind; 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 larger, alternately changing the direction of the streamlines; at T=15000, two large vortices are formed on the left and right sides, and the direction of the streamlines no longer changes with the wind. These simulation results are consistent with real physical laws, proving the effectiveness of the method in simulating large-scale hydrodynamic models of the Earth's surface.

[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 velocity and inundation area of ​​the flood at different times were analyzed. A schematic diagram of the bottom slope height of the river section in the dynamic simulation model of the dam-break flood is shown below. Figure 7 As shown.

[0104] By solving the Boltzmann equations in shallow water lattice and importing measured data from river cross-sections, an adaptive grid technique is employed to achieve a refined representation of topographic features. In a hypothetical dam failure scenario at a cascade hydropower station downstream of the Zipingpu Reservoir, a breach width of 200 meters and an initial peak flow rate of 9000 m³ / h are set. 3 / s, to perform dynamic simulation of the flood evolution process for up to 1 hour. For example... Figure 8 As shown, this plot displays flood flow velocity maps at different times, including the inundated areas.

[0105] Simulation results show that the dam-break wave generates an impact velocity of up to 5 m / s in the narrow river section, causing some low-lying areas to be flooded within 40 minutes after the dam break, with a maximum water depth of 4.5 meters. By outputting parameters such as flood depth, flow velocity vector field, and flood arrival time, the interaction mechanism between the dam-break flood and the terrain is analyzed.

[0106] The method of this invention took 23 minutes to run on a computer equipped with an Intel i7-12700 processor, while it took 51 minutes on the same platform using Mike 21, a well-known industry-leading flood simulation software. Engineering applications show that the method of this invention not only significantly reduces computation time in flood evolution simulation in complex terrain, but also improves the predictive reliability of key parameters such as inundation 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 hydrodynamic simulation method based on shallow water equations, which significantly improves the computational stability of dynamic boundary scenarios such as dam break wavefront and intertidal zone inundation while ensuring mass-momentum conservation, and shortens the computation time.

[0108] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions 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 large-scale surface hydrodynamic simulation based on shallow water equations, characterized by the following steps: include: S1. Introduce the bottom source term to construct the shallow water lattice Boltzmann equation; S2. Constructing equations for handling wet and dry boundaries using a single relaxation model, specifically including: By combining a single relaxation model to construct equations for treating wet and dry boundaries, in the shallow water lattice Boltzmann equations, the flow passes through the distribution function. The evolution of collision and migration processes is expressed by the following formula: ； In the formula, External force term; The distribution function is decomposed into equilibrium and non-equilibrium disturbance terms, as shown in the formula: ； In the formula, ϵ is a dimensionless small parameter; at the dry element... h =0, velocity u =0, This is a non-equilibrium term, propagated through spatial gradient; Discretize the spacetime, perform Taylor expansion on time and space, and combine the discretization scheme to approximate the partial derivatives. The formula is as follows: ； In the formula, Let h represent the distribution function of the wet element, where h > 0 at the wet element. Represents the distribution function of the dry unit; Calculate the distribution function of the wetted element The spatial gradient, based on the wet unit distribution function Spatial ladder calculation dry element distribution function ,include: Calculate the distribution function of the wetted element The formula for the spatial gradient is: ； In the formula, For discrete velocity direction index; For spatial step size, The distribution function value is the value on the right side of the right wetted element. The value of the distribution function on the left side of the wet cell; Dry unit distribution function The formula after spatial gradient correction is: ； In the formula, For correction factor, For gradient information; S3. Construct the model to be simulated and set the computational domain according to the constructed model. The model to be simulated is a dam break wave scour cone group and slope revetment model, including upstream, downstream and downstream slope regions. The downstream region includes the cone group. The computational domain includes two variables: length and width. S4. Set the simulation parameters of the model to be simulated, including the calculation time step, fluid density, relaxation factor, gravitational acceleration, Manning coefficient and wind speed; S5. Set the boundary conditions for the model to be simulated. Set the upstream region and the left and right sides of the flow field as non-slip solid wall boundaries, and set the cone group and the downstream slope region as non-equilibrium extrapolation boundary conditions. S6. Simulate the dynamic phenomena of the model to be simulated based on the constructed shallow water lattice Boltzmann equation and the wet and dry boundary treatment equation.

2. The large-scale surface hydrodynamic simulation method based on shallow water equations according to claim 1, characterized in that, Step S1 specifically includes: The Boltzmann equations for shallow-water lattice systems are constructed as follows: ； In the formula, It is the particle distribution function; It is a spatial vector, in two-dimensional space. ; It is time; It is the time step; It is the particle velocity vector; yes In the Components in each direction; , It's the size of the grid; It is the elevation of the bed surface; It is a single relaxation time; It is a local equilibrium distribution function; Indicates the weighting coefficient; ； ； In the formula, It is the volume force acting on the fluid. directional components, It refers to the water depth, which is the height from the bed surface to the water surface. It is the depth-average velocity; It is the first Wind shear stress in each direction; It is the shear stress on the bed surface; It is the density of water; The depth-averaged velocity is at the 1st Components in direction; Is the depth-average velocity at y Components in direction; Is the depth-average velocity at x Components in direction; It is the coefficient of friction of the bed surface; It is the Coriolis parameter; It is the Kronecker function, and the formula is: ； The general local equilibrium distribution function is redefined as follows: ； In the formula, A and B It is an unknown constant; g It is gravitational acceleration; It is the particle velocity vector e α exist i The directional component.

Images