A simulation calculation method for two-dimensional windblown sand transport rate on the coastal dune plane

Abstract

The present invention relates to a simulation calculation method for two-dimensional wind-blown sand transport rate on the coastal dune plane. It reads the coastal terrain, boundaries, planar triangular meshes, and keyword files to obtain coastal terrain data, representative sediment particle sizes, and wind speed data at the boundaries; continues to obtain the critical starting wind speed, derives the elevation gradient, and performs interpolation calculations; calculates the wind-blown sand transport rate. The present invention realizes the calculation of wind speed changes, wind direction deflection, and transport rate generated when the wind incident on the beach at a certain angle and a certain height passes through terrains with different elevations on the beach surface.

Description

Technical Field

[0001] The present invention relates to a simulation calculation method for the planar two-dimensional wind-sand sediment transport rate of coastal dunes, belonging to the fields of coastal dune protection and restoration and wind-sand prevention and control. Background Art

[0002] With global climate change and sea-level rise, the risks of overwash and flood caused by storm surges have increased, leading to a sudden increase in the pressure of coastal disaster prevention and mitigation. Coastal dunes are the last line of defense against coastal floods and play an important role in protecting the lives and property of people in coastal areas. The protection and restoration of sandy coastal dunes have become one of the most important construction contents of the coastal zone protection and restoration project in China.

[0003] Coastal wind-sand movement refers to the movement of near-bottom sediment mainly in a saltating manner driven by the wind, which is the main process of sediment exchange in coastal dunes. Under different dynamic conditions, wind-sand movement causes rapid (or slow) changes in the topography of coastal dunes. Although during storm surges, the sediment movement along the coast is mainly water-sediment movement, during normal waves (or storm recovery intervals), wind-sand movement plays a dominant role. Accurate calculation of the wind-sand sediment transport rate is a prerequisite for successfully simulating the topographic evolution of coastal dunes and an inevitable requirement for the protection and restoration of coastal dunes. However, coastal wind-sand movement is complex. On the time scale, due to the high randomness of wind power, wind-sand movement is highly dynamic; on the spatial scale, due to the complex topography of dunes and rich geomorphic elements on the actual coast, the coupling effect with the wind field makes the calculation of the wind-sand sediment transport rate more cumbersome.

[0004] The bottleneck in the simulation of wind-sand movement at the present stage is that most of the current experiments are still focused on simulating the wind-sand sediment transport rate along a specific cross-section direction under the action of the incident wind perpendicular to the shoreline [1] , while the wind on the actual coast often incides at a certain angle with the shoreline, which involves the process of coastal sediment exchange, and there is no corresponding simulation calculation method for the wind-sand sediment transport rate in a planar area. For the wind incident at a certain angle with the shoreline on the actual coast, the existing domestic and foreign models only simulate the two-dimensional wind field incident on the dunes at a certain angle, without involving the calculation of the sediment transport rate or only simulating the sediment transport rate at different positions under a one-dimensional wind field. The present invention mainly aims at the above problems and proposes a simulation calculation method for the planar two-dimensional wind-sand sediment transport rate of coastal dunes.

[0005] [1]Roelvink,D.,&Costas,S.(2019).Coupling nearshore and aeolian processes:XBeach and duna process-based models.Environmental Modelling&Software,115,98-112. Summary of the Invention

[0006] The present invention provides a simulation calculation method for the two-dimensional sand-drift rate on the plane of coastal dunes, which realizes the calculation of the wind speed change, wind direction deflection and sand-drift rate generated when the wind incident on the beach at a certain angle and height passes through the terrain at different elevations on the beach surface.

[0007] The technical solution adopted by the present invention to solve its technical problems is as follows:

[0008] A simulation calculation method for the two-dimensional sand-drift rate on the plane of coastal dunes specifically includes the following steps:

[0009] Step S1: Read the coastal terrain, boundaries, planar triangular meshes and keyword files to obtain coastal terrain data, representative particle sizes of sediment, and wind speed data at the boundaries;

[0010] Step S2: Set the x-direction to represent the offshore direction, and the offshore direction is the positive direction; the y-direction represents the longshore direction, and the downstream direction is the positive direction; the z-direction is the terrain elevation direction, and vertically upward is the positive direction; assume that the wind at the boundary is incident on the coast at a preset angle incident on the coast, where is the angle formed by the incident direction and the positive x-axis direction;

[0011] The formula representing the distribution of wind shear stress on the terrain profile is:

[0012]

[0013] In Formula (1) and Formula (2), represents the surface shear stress perturbation value in the x-direction, represents the surface shear stress perturbation value in the y-direction, represents the partial derivative of the ground elevation in the x-direction, represents the partial derivative of the ground elevation in the y-direction;

[0014] Step S3: Calculate the terrain elevation gradient using planar triangular meshes and

[0015] Step S4: Based on the terrain elevation gradient obtained in Step S3, make profiles along the x-axis and y-axis respectively at the planar triangular grid nodes (x, y) where the wind speed is to be calculated. Take equally spaced interpolation nodes on the profiles, and interpolate and respectively at the interpolation nodes to obtain the terrain elevation gradient at any position of the planar triangular grid; and respectively, so as to obtain the terrain elevation gradient at any position of the planar triangular grid;

[0016] Step S5: Denote any interpolation node on the profile as (x0, y0), select the four grid nodes (x i , y i ) closest to (x0, y0), and calculate their average value as the function value at the interpolation node;

[0017] Step S6: Based on the function values at the interpolation nodes obtained in Step S5, integrate formulas (1) and (2) on the x-profile and y-profile to calculate the wind speed components at a certain height of the point to be determined (x, y):

[0018]

[0019] Here, the certain height of the point to be determined (x, y) refers to the height equal to that of the anemometer for measuring the wind speed on the actual coast;

[0020] Step S7: Calculate the near-bottom critical incipient velocity according to the representative sediment grain size obtained in Step S1:

[0021]

[0022] In formula (10), u* t is the near-bottom critical incipient velocity, A is an empirical coefficient (take 0.08), ρ p [kg / m 3 is the sediment particle density, ρ a [kg / m 3 is the air density, d n [m] is the representative sediment grain size;

[0023] Step S8: Convert the near-bottom critical incipient velocity to the critical incipient velocity at a certain height through the wall function:

[0024] u t = bu *t (11)

[0025]

[0026] In formulas (11) and (12), z’ is the near-bottom roughness height, that is, the height where the theoretical wind speed is 0, and κ is the von Kármán constant, take 0.4;

[0027] Step S9: Calculate the wind-sand sediment transport rate:

[0028]

[0029] where C is the particle size distribution width parameter, g [m / s 2 is the acceleration due to gravity, D [m] is the reference particle size; Q [kg / m / s] is the sediment transport rate, u [m / s] is the wind speed at height z, and u t [m / s] is the critical starting wind speed at the same height;

[0030] As a further preference of the present invention,

[0031] In step S2, α and β are coefficients related to the terrain elevation characteristic parameters and the near-bottom roughness height, taking 3 and 0.2 respectively for the windward slope and 2 and 1 respectively for the leeward slope;

[0032] As a further preference of the present invention,

[0033] In step S3, the planar triangular grid element has three vertices, represented by node numbers i = 1, 2, 3 or j = 1, 2, 3. Then, in the planar triangular grid element, the elevation difference between adjacent nodes is defined as follows:

[0034]

[0035] In formula (3), the subscripts i and j represent the coordinate value of element number i minus the coordinate value of element number j;

[0036] The elevation gradient of the element nodes and The calculation formula is:

[0037]

[0038] In formula (4), the subscripts i, j, and k are the element node numbers; they are taken in positive order in sequence, for example, i = 1, j = 2, k = 3 or i = 2, j = 3, k = 1;

[0039] As a further preference of the present invention,

[0040] In step S5, the interpolation method is specifically: Denote any interpolation node on the section as (x0, y0), and select the four grid nodes (x i , y i ) closest to (x0, y0), where i = 1, 2, 3, 4, and the four selected grid nodes satisfy:

[0041]

[0042] And for the and Calculate the average values respectively as the function values at the interpolation nodes, and obtain:

[0043]

[0044] In formula (6), and are the values at the interpolation nodes.

[0045] By the above technical solutions, compared with the prior art, the present invention has the following beneficial effects:

[0046] The simulation calculation method for the planar two-dimensional wind-sand sediment transport rate of coastal dunes provided by the present invention can simulate the influence of wind direction deflection on the spatial distribution of the sediment transport rate under two-dimensional terrain, and realize the calculation of the wind speed change, wind direction deflection and sediment transport rate generated when the wind incident on the beach at a certain angle and a certain height passes through terrains with different elevations on the beach surface. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] The present invention will be further described below in conjunction with the drawings and embodiments.

[0048] Figure 1 is a schematic flow chart of the simulation calculation method for the planar two-dimensional wind-sand sediment transport rate of coastal dunes provided by the present invention;

[0049] Figure 2 is the change and deflection of the wind speed when the wind at different angles with the positive direction of the x-axis is incident on the dune from the position of x = -83m in the preferred embodiment provided by the present invention. Among them, the direction indicated by the red arrow represents the wind direction, and the depth of the topographic map color represents the magnitude of the wind speed;

[0050] Figure 3 is a representative longitudinal profile of the wind speed change selected in the preferred embodiment provided by the present invention, showing the change trend of the wind speed at different incident angles affected by the ground elevation;

[0051] Figure 4 is the spatial distribution of the wind-sand sediment transport rate caused by the wind at different incident angles in the preferred embodiment provided by the present invention;

[0052] Figure 5 is a comparison diagram of the wind speed simulated on a typical one-dimensional profile by the method proposed in the present invention and the wind speed obtained by running the same example on the XBeach Duna model in the preferred embodiment provided by the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0053] The present invention will now be described in further detail with reference to the accompanying drawings. In the description of this application, it should be understood that the orientation or positional relationships indicated by terms such as "left side", "right side", "upper part", "lower part", etc. are based on the orientation or positional relationships shown in the drawings. These are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the devices or elements referred to must have a specific orientation, be constructed and operated in a specific orientation. "First", "second", etc. do not represent the importance of components, so they should not be construed as limitations on the present invention. The specific dimensions used in this embodiment are only for illustrative purposes of the technical solution and do not limit the protection scope of the present invention.

[0054] Coastal dunes are the last line of defense against coastal floods and play an important role in protecting the lives and property of people in coastal areas. Accurately calculating the sand transport rate is a prerequisite for successfully simulating the topographic evolution of coastal dunes and an inevitable requirement for the protection and restoration of coastal dunes. However, coastal aeolian sand movement is complex. On the time scale, due to the high randomness of wind power, the aeolian sand movement has a high degree of dynamics; on the spatial scale, due to the complex dune topography and rich geomorphic elements on the actual coast, the coupling effect with the wind field makes the calculation of the sand transport rate more cumbersome.

[0055] As described in the background art, at present, the simulation of aeolian sand movement mainly focuses on simulating the sand transport rate along a specific profile direction under the action of incident wind perpendicular to the shoreline. However, the wind on the actual coast often incides at a certain angle with the shoreline, which involves the process of coastal sediment exchange, and there is no corresponding simulation calculation method for the sand transport rate in a planar area. Even though there are a few models at home and abroad that can simulate the two-dimensional wind field incident on the dunes at a certain angle, they do not involve the calculation of the sand transport rate or only simulate the sand transport rate at different positions under a one-dimensional wind field. Therefore, this application attempts to overcome the problems that cannot be considered by the existing one-dimensional sand transport rate calculation models, that is, to simulate the influence of wind direction deflection on the spatial distribution of the sand transport rate under two-dimensional topography.

[0056] As Figure 1 shown, it is a simulation calculation method for the planar two-dimensional sand transport rate of coastal dunes provided by the present invention, specifically including the following steps:

[0057] Step S1: Read the coastal topography, boundary, planar triangular grid, and keyword file to obtain the coastal topography data, representative particle size of sediment, and wind speed data at the boundary.

[0058] Step S2: Set the x-direction to represent the direction away from the shore, and the direction away from the shore is the positive direction; the y-direction represents the longshore direction, and the downstream direction is the positive direction; the z-direction is the terrain elevation direction, and vertically upward is the positive direction; assume that the wind at the boundary incides on the coast at a preset angle incides on the coast, is the angle formed by the incident direction and the positive direction of the x-axis, and

[0059] The formula for characterizing the distribution of wind shear stress on the terrain profile is:

[0060]

[0061] In Formula (1) and Formula (2), represents the surface shear stress perturbation value in the x direction, represents the surface shear stress perturbation value in the y direction, represents the partial derivative of the ground elevation in the x direction, represents the partial derivative of the ground elevation in the y direction. α and β are coefficients related to the terrain elevation characteristic parameters and the near-bottom roughness height, taking 3 and 0.2 respectively for the windward slope and 2 and 1 respectively for the leeward slope; the above two formulas are based on the incompressible turbulent boundary layer flow and are derived to simulate the wind speed change formula of the turbulent wind passing over the sand dune.

[0062] Step S3: Calculate the terrain elevation gradient using a planar triangular grid and A planar triangular grid element has three vertices, represented by node numbers i = 1, 2, 3 or j = 1, 2, 3. Then, in the planar triangular grid element, the elevation difference between adjacent nodes is defined as follows:

[0063]

[0064] In Formula (3), the subscripts i and j represent the coordinate value of element number i minus the coordinate value of element number j;

[0065] The elevation gradient of the element nodes and The calculation formula is:

[0066]

[0067] In Formula (4), the subscripts i, j, and k are the element node numbers; they are taken in positive order in sequence, for example, i = 1, j = 2, k = 3 or i = 2, j = 3, k = 1.

[0068] This step reflects the innovation degree of the invention point of this application. By starting from calculating the gradients of the terrain elevation of each node in the x direction and the y direction, the purpose of calculating the entire two-dimensional plane is achieved.

[0069] Step S4: Based on the terrain elevation gradient obtained in Step S3, make profiles along the x-axis and the y-axis respectively at the planar triangular grid nodes (x, y) where the wind speed is to be calculated, take equally spaced interpolation nodes on the profiles, and at the interpolation nodes for Interpolate separately with to obtain the terrain elevation gradient at any position of the planar triangular mesh.

[0070] Step S5: Denote any interpolation node on the cross-section as (x0, y0), and select the four grid nodes (x i , y i ) closest to (x0, y0), where i = 1, 2, 3, 4, and the four selected grid nodes satisfy:

[0071]

[0072] And calculate the average values of and on the 4 grid nodes respectively as the function values at the interpolation node, obtaining:

[0073]

[0074] In formula (6), and are the values at the interpolation node.

[0075] The main purpose of this step is to derive the partial derivatives of the ground elevation in the x-direction and y-direction at the cross-section interpolation node, so as to be able to integrate formulas (1) and (2) on the x-cross-section and y-cross-section to calculate the wind speed components at the point to be determined.

[0076] Step S6: Based on the function values at the interpolation nodes obtained in step S5, integrate formulas (1) and (2) on the x-cross-section and y-cross-section to calculate the wind speed components at a certain height of the point to be determined (x, y):

[0077]

[0078] The certain height of the point to be determined (x, y) here refers to the height equal to that of the anemometer measuring the wind speed on the actual coast.

[0079] Here, a specific elaboration is needed. The reason for calculating the wind speed components is that the change in the wind speed components will cause changes in the total wind speed and total wind direction. Formula (9) is the change formula for the total wind speed, which is determined by the wind speed components in the x-direction and y-direction and changes with the numerical values of the wind speed components in the x-direction and y-direction. The wind obliquely incident on the dune from the incident boundary. During the process of the wind passing over the dune, due to the spatial change of the ground elevation, the wind speed components also change spatially (this conclusion can be obtained from formulas (1), (2), (7), and (8)), resulting in a change in the total wind direction and gradually becoming perpendicular to the dune crest line. This is the phenomenon of wind deflection. Therefore, calculating the wind speed components is necessary for calculating the total wind speed, total wind direction, and wind deflection.

[0080] Step S7: Calculate the near-bottom critical incipient velocity based on the representative sediment particle size obtained in Step S1:

[0081]

[0082] In Equation (10), u* t is the near-bottom critical incipient velocity, A is an empirical coefficient (take 0.08), ρ p [kg / m 3 is the sediment particle density, ρ a [kg / m 3 is the air density, d n [m] is the representative sediment particle size.

[0083] Step S8: Convert the near-bottom critical incipient velocity to the critical incipient velocity at a certain height through the wall function:

[0084] u t = bu *t (11)

[0085]

[0086] In Equations (11) and (12), z’ is the near-bottom roughness height, i.e., the height where the theoretical wind speed is 0, and κ is the von Kármán constant, take 0.4.

[0087] Step S9: If the critical incipient velocity of the point to be determined (x, y) is greater than the wind speed of the point to be determined, the sand-dust transport rate shown in the output result file is 0. If the critical incipient velocity of the point to be determined (x, y) is less than the wind speed of the point to be determined, calculate the sand-dust transport rate using the following formula:

[0088]

[0089] Among them, C is the particle size distribution width parameter, g [m / s 2 is the gravitational acceleration, D [m] is the reference particle size; Q [kg / m / s] is the sediment transport rate, u [m / s] is the wind speed at height z, and u t [m / s] is the critical incipient velocity at the same height.

[0090] To verify the excellence of this application, this application provides an example of deriving the sand-dust transport rate in a two-dimensional plane using the above method, and compares its calculation results with those of the one-dimensional profile model (XBeach Duna, reference 1) of the sand-dust transport rate commonly used at present. The simulated terrain is a certain two-dimensional sand dune, the incident wind direction forms an angle of 30 - 60 degrees with the positive x-axis, the incident wind speed is defined as the wind speed at a height of 10 m above the ground, with a magnitude of 10 m / s, and it is incident from the position of x = -83 m.

[0091] Figure 2It shows the change and deflection of the wind speed when the wind blows into the sand dune from the position of x = -83 m at different angles with the positive direction of the x-axis. The direction indicated by the red arrow represents the wind direction, and the depth of the topographic map color represents the magnitude of the wind speed.

[0092] Figure 3 Typical longitudinal profiles of wind speed changes are selected to more clearly describe the changing trends of wind speeds at different incident angles under the influence of ground elevation.

[0093] Figure 4 It shows the spatial distribution of the sand transport rate caused by winds at different incident angles.

[0094] Figure 5 This is a comparison chart of the wind speed distribution along a one-dimensional characteristic profile arbitrarily selected along the x-axis for simulating the wind speed distribution under a two-dimensional terrain using the method proposed in the present invention and the wind speed distribution of the one-dimensional characteristic profile simulated using the XBeach Duna model under the same incident wind conditions. It can be seen that the error between the wind speed distribution of any one-dimensional profile under the two-dimensional terrain simulated using the present invention and the wind speed distribution of the same one-dimensional profile simulated by XBeach Duna is within a reasonable range, and it can correctly reflect the changing trend of the wind speed, indicating that the accuracy and reliability of the simulation of the two-dimensional wind speed distribution by the present invention can meet the actual requirements.

[0095] Those skilled in the art of this technology can understand that, unless otherwise defined, all terms (including technical terms and scientific terms) used here have the same meaning as the general understanding of those of ordinary skill in the technical field to which this application belongs. It should also be understood that terms such as those defined in a general dictionary should be understood to have a meaning consistent with the meaning in the context of the prior art, and will not be interpreted with an idealized or overly formal meaning unless defined as here.

[0096] The meaning of "and / or" described in this application refers to the situation where each exists alone or both exist simultaneously.

[0097] The meaning of "connection" described in this application can be a direct connection between components or an indirect connection between components through other components.

[0098] Taking the ideal embodiments of the present invention as the inspiration, through the above description, relevant staff can completely make various changes and modifications without departing from the technical idea of this invention. The technical scope of this invention is not limited to the content in the specification, and its technical scope must be determined according to the scope of the claims.

Claims

1. A simulation calculation method for two-dimensional wind-blown sand transport rate on the coastal dune plane, characterized in that: Specifically, it includes the following steps: Step S1: Read the coastal terrain, boundary, planar triangular mesh, and keyword files to obtain coastal terrain data, representative sediment particle sizes, and wind speed data at the boundary; Step S2: Set the x-direction to represent the direction towards the offshore, with the offshore direction as the positive direction; the y-direction to represent the longshore direction, with the downstream direction as the positive direction; the z-direction to be the terrain elevation direction, with vertically upward as the positive; assume that the wind at the boundary is incident on the coast at a preset angle onto the coast, where the angle is the angle formed by the incident direction and the positive x-axis direction; The formula for characterizing the distribution of wind shear stress on the terrain profile is as follows: In Formula (1) and Formula (2), represents the surface shear stress perturbation value in the x direction, represents the surface shear stress perturbation value in the y direction, represents the partial derivative of the ground elevation in the x direction, represents the partial derivative of the ground elevation in the y direction; Step S3: Calculate the terrain elevation gradient using a planar triangular grid and Step S4: Based on the terrain elevation gradient obtained in Step S3, make profiles along the x-axis and y-axis respectively at the planar triangular grid nodes (x, y) where the wind speed is to be calculated. Take equally spaced interpolation nodes on the profiles, and interpolate and respectively to obtain the terrain elevation gradient at any position of the planar triangular grid; Step S5: Denote any interpolation node on the section as (x0, y0), and select the four grid nodes (x i , y i ) closest to (x0, y0), and calculate their average value as the function value at the interpolation node; Step S6: Based on the function values at the interpolation nodes obtained in Step S5, integrate formulas (1) and (2) on the x-section and y-section to calculate the wind speed components at a certain height at the point to be determined (x, y): The certain height at the point to be determined (x, y) here refers to the height equal to that of the anemometer for measuring wind speed on the actual coast; Step S7: Calculate the near-bottom critical starting wind speed according to the representative sediment particle sizes obtained in Step S1; In formula (10), u* t is the near-bottom critical starting wind speed, A is an empirical coefficient, taking 0.08, ρ p [kg / m 3 is the sediment particle density, ρ a [kg / m 3 is the air density, d n [m] is the representative sediment particle size; Step S8: Convert the near-bottom critical starting wind speed to the critical starting wind speed at a certain height through the wall function; u t = bu *t (11) In formulas (11) and (12), z’ is the near-bottom roughness height, that is, the height where the theoretical wind speed is 0, and κ is the von Kármán constant, taking 0.4; Step S9: Calculate the wind-blown sand transport rate; where C is the particle size distribution width parameter, g [m / s 2 is the acceleration due to gravity, D [m] is the reference particle size; Q sat [kg / m / s] is the sediment transport rate, u [m / s] is the wind speed at height z, u t [m / s] is the critical incipient motion wind speed at the same height.

2. The simulation calculation method of the planar two-dimensional wind-blown sand transport rate of coastal dunes according to claim 1, wherein: In Step S2, α and β are coefficients related to the terrain elevation characteristic parameters and the near-bottom roughness height, taking 3 and 0.2 respectively for the windward slope, and 2 and 1 respectively for the leeward slope.

3. The simulation calculation method of the planar two-dimensional wind-blown sand transport rate of coastal dunes according to claim 1, wherein: In Step S3, the planar triangular mesh element has three vertices, represented by node numbers i = 1, 2, 3 or j = 1, 2, 3. Then, in the planar triangular mesh element, the elevation differences between adjacent nodes are defined as follows: In formula (3), the subscripts i and j represent the coordinate value of element number i minus the coordinate value of element number j; Elevation gradient of unit nodes and The calculation formula is as follows: In formula (4), the subscripts i, j, and k are the element node numbers; taking values in positive order in turn, i = 1, j = 2, k = 3 or i = 2, j = 3, k = 1.

4. The simulation calculation method of the planar two-dimensional wind-blown sand transport rate of coastal dunes according to claim 1, wherein: In step S5, the interpolation method is specifically as follows: Denote any interpolation node on the section as (x0, y0), and select the four grid nodes (x i , y i ) that are closest to (x0, y0), where i = 1, 2, 3, 4, and the four selected grid nodes satisfy: and calculate the average values for the and at the four grid nodes respectively as the function values at the interpolation nodes, obtaining: In formula (6), and are the values at the interpolation nodes.

Images