CN109582996A - A kind of coupled simulation method of small scale beach profile and large scale Shoreline changes - Google Patents
A kind of coupled simulation method of small scale beach profile and large scale Shoreline changes Download PDFInfo
- Publication number
- CN109582996A CN109582996A CN201810944700.4A CN201810944700A CN109582996A CN 109582996 A CN109582996 A CN 109582996A CN 201810944700 A CN201810944700 A CN 201810944700A CN 109582996 A CN109582996 A CN 109582996A
- Authority
- CN
- China
- Prior art keywords
- wave
- formula
- scale
- water
- model
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 65
- 238000004088 simulation Methods 0.000 title claims abstract description 27
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims abstract description 83
- 239000013049 sediment Substances 0.000 claims abstract description 58
- 239000004576 sand Substances 0.000 claims description 62
- 230000008859 change Effects 0.000 claims description 45
- 230000008569 process Effects 0.000 claims description 33
- 238000004364 calculation method Methods 0.000 claims description 31
- 230000004907 flux Effects 0.000 claims description 27
- 238000009826 distribution Methods 0.000 claims description 25
- 230000008878 coupling Effects 0.000 claims description 24
- 238000010168 coupling process Methods 0.000 claims description 24
- 238000005859 coupling reaction Methods 0.000 claims description 24
- 238000009792 diffusion process Methods 0.000 claims description 18
- 239000000126 substance Substances 0.000 claims description 15
- 230000033001 locomotion Effects 0.000 claims description 14
- 238000012937 correction Methods 0.000 claims description 12
- 230000003628 erosive effect Effects 0.000 claims description 10
- 239000002245 particle Substances 0.000 claims description 10
- 230000005855 radiation Effects 0.000 claims description 8
- 239000000463 material Substances 0.000 claims description 7
- 239000002352 surface water Substances 0.000 claims description 5
- 230000008021 deposition Effects 0.000 claims description 4
- 238000004062 sedimentation Methods 0.000 claims description 4
- 230000001133 acceleration Effects 0.000 claims description 3
- 239000012530 fluid Substances 0.000 claims description 3
- 238000011010 flushing procedure Methods 0.000 claims description 3
- 238000013467 fragmentation Methods 0.000 claims description 3
- 238000006062 fragmentation reaction Methods 0.000 claims description 3
- 230000005484 gravity Effects 0.000 claims description 3
- 230000008595 infiltration Effects 0.000 claims description 3
- 238000001764 infiltration Methods 0.000 claims description 3
- 230000010355 oscillation Effects 0.000 claims description 3
- 230000021715 photosynthesis, light harvesting Effects 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 238000010008 shearing Methods 0.000 claims 2
- 230000001419 dependent effect Effects 0.000 claims 1
- 230000001808 coupling effect Effects 0.000 abstract description 5
- 238000011002 quantification Methods 0.000 abstract 2
- 230000007774 longterm Effects 0.000 description 10
- 230000009471 action Effects 0.000 description 9
- 230000003993 interaction Effects 0.000 description 7
- 230000000694 effects Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 238000009991 scouring Methods 0.000 description 5
- 238000005096 rolling process Methods 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 3
- 230000009467 reduction Effects 0.000 description 3
- 238000012795 verification Methods 0.000 description 3
- 238000006424 Flood reaction Methods 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 2
- 230000003247 decreasing effect Effects 0.000 description 2
- 239000006185 dispersion Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000000746 purification Methods 0.000 description 2
- 208000011117 substance-related disease Diseases 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 239000000725 suspension Substances 0.000 description 2
- 238000009825 accumulation Methods 0.000 description 1
- 230000001154 acute effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000010924 continuous production Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000005137 deposition process Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000013178 mathematical model Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000002265 prevention Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000000630 rising effect Effects 0.000 description 1
- 229920006395 saturated elastomer Polymers 0.000 description 1
- 230000007480 spreading Effects 0.000 description 1
- 238000003892 spreading Methods 0.000 description 1
- 239000004575 stone Substances 0.000 description 1
- 238000001363 water suppression through gradient tailored excitation Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses the approach of coupled numerical simulation of a kind of small scale beach profile and large scale Shoreline changes, include the following steps: the Coupling effect model that S1, foundation describe plane offshore power accurate quantification；The bottom offshore flow model that S2, foundation describe offshore bottom offshore current accurate quantification；S3, small scale beach profile evolution model is established；S4, using long period control function sets, establish large scale Shoreline changes model；S5, using Godunov type finite volume method on unstrctured grid, solve large scale Shoreline changes model, obtain large scale Shoreline changes numerical value.The present invention is realized from the expansion of Sediment Transport amount conveying capacity on large scale water front landforms on small scale seashore, is completed the simulation of large scale Shoreline changes, is remarkably improved the practical application value of beach evolution analog simulation.
Description
Technical Field
The invention relates to a numerical simulation method for beach evolution of a sandy and silty estuary coast, belonging to the technical field of numerical simulation of movement of silt and beach evolution of the estuary coast.
Background
The coastal region is a region where the economy and the society of China develop at a high speed, and has a strategic position of putting a great importance in the economic construction of China. However, the influence of people on the landform of coastal areas is increasingly intensified, and particularly, the sustainable utilization and protection of coastal zone resources, shipping industry, flood prevention and disaster reduction, ecological environment protection and other important problems are influenced for a long time. The method is limited by field conditions (flood, storm surge), observation instruments, capital and the like, and the method becomes an important technical means for research by utilizing a mathematical model facing the coastal dynamic landform process to predict the long-term growth-diminishing and short-term evolution law of the coastline.
The river mouth and coast landform evolution is an extremely complex process, comprises the coast short-term adjustment and medium-and-long-term erosion evolution processes, and has a multi-space-time scale phenomenon. The physical process influencing the evolution of the landform of the estuary and the coast not only has direct factors such as water flow, wave and sediment transportation, but also comprises long-term factors such as incoming water, incoming sand, wind sand transportation and the like, and the factors occur at different time and space scales, have nonlinear influence on the coast form and need to be organically combined. Thus, model building requires efficient coupling of functional modules that describe physical processes at different scales and computing them at scales corresponding to the physical processes. In order to simulate the evolution of the coastal system more truly, the interaction of different space-time scales and landform units of the coastal system needs to be further coupled, and the problems of the evolution and scale matching are solved.
In the traditional estuary and coastal landform evolution simulation, a specific development process of a small-scale phenomenon is generally researched in the face of process simulation, and a real process is ignored in the face of performance simulation, and only the whole change of a variable scale is concerned. Therefore, the model with the coupling effect of power factors with different space-time scales neglected has poor universality, and is difficult to meet the actual requirements of sustainable utilization and protection of coastal zones, flood control and disaster reduction, ecological environment protection and the like.
Disclosure of Invention
The invention aims to overcome the defects and shortcomings of the prior art, provides a numerical simulation method for coupling small-scale beach section and large-scale shoreline change based on a sand conveying process, realizes the description on the coupling of a estuary coast landform evolution dynamic process, solves the problems of matching of multiple space-time scales in evolution, calculating time consumption in long-period simulation and the like, and obviously improves the practical application value of beach evolution simulation.
In order to achieve the purpose of the invention, the invention adopts the following technical scheme: a coupling simulation method for small-scale beach profile and large-scale shoreline change comprises the following steps:
s1, establishing a wave-current coupling model for calculating main hydrodynamic factors (such as tide, wind-driven circulation, runoff, near-shore current caused by wave breaking and the like) influencing the beach evolution so as to realize accurate quantitative description of the plane near-shore power;
s2, establishing a bottom offshore flow model; the step takes the wave element provided by S1 as the input condition, the turbulent fluctuation intensity is expressed based on the linear vortex viscosity coefficient vertical distribution, the bottom off-shore flow is driven by the wave residual momentum flow gradient and the pressure gradient caused by increasing and decreasing the water level together, a bottom off-shore flow model is established for calculating the bed shear stress and the water flow velocity distribution at any time under the coupling action of wave and flow, the accurate quantitative tracing of the off-shore flow at the near-shore bottom is realized,
s3, establishing a small-scale beach profile evolution model for calculating sediment transport under the action of near shore power and beach change caused by sediment movement, and realizing accurate quantitative description of interaction of a hydrodynamic process, a sediment transport process and beach landform evolution;
s4, establishing a large-scale shoreline change model by using the long-period control function set, and expanding the scale in the small-scale shoal section model to realize quantitative prediction of large-scale shoreline change;
and S5, solving the large-scale shoreline change model by adopting a Godunov finite volume method on the unstructured grid, and acquiring a large-scale shoreline change value.
Further, the step S1 specifically includes the following steps:
s11, establishing a two-dimensional hydrodynamic module: in the invention, a conservation-form two-dimensional shallow water equation is used as a water flow control equation, which is as follows:
in the formula,is a conservation vector;、are respectively as、A directional convective flux vector;、are respectively as、Diffusion flux vector caused by directional reynolds stress;、are respectively as、A diffusion flux vector caused by directional secondary flow;for the source term vector:
；；；
；；
；；
；
in the formula,is the water depth;、respectively, average flow velocity in the vertical direction、A component of direction;is the bottom elevation;the intensity of rainfall is used;the infiltration strength is obtained;is the turbulent viscosity coefficient in the horizontal direction,，is a coefficient of proportionality that is,in order to be the karman coefficient,the flow rate is cut for the bed surface;、、、is a diffusion stress term caused by secondary flow;is the acceleration of gravity;is coefficient of Coriolis force，，Is the angular velocity of the earth's rotation,the local latitude is;in order to be the wind stress,、，the water surface wind stress drag coefficient;、density of air and water, respectively;、are respectively as、A wind speed component (m/s) at a height of 10m above the water surface in the direction;、、the wave radiation stress is provided by a wave model;、are respectively as、The slope of the friction resistance in the direction,、. Wherein,h、u、vis the amount of strain.
S12, establishing a wave module, wherein the wave module adopts the third generation wave of SWAN based on the conservation equation of wave action quantity. A wave field module for providing wave elementsE w 、E r 、、HAnd calculating wave radiation stress and wave surface water rolling momentum flow. The specific formula is as follows:
in the formula:is the wave propagation direction andthe included angle of the axes;is the average wave energy per wave period of the water column,，nis the ratio of the group velocity to the phase velocity,，in terms of the wave number, the number of waves,E r is the energy density of the water roll.
Because the wave radiation stress changes along the way, the mean water surface changes, and the mean water surface elevation (water increase and decrease) can be calculated by the following formula in combination with the consideration of the influence of the rolling water:
and S13, coupling the two-dimensional hydrodynamic module in the S11 and the wave module in the S12 to form a wave current coupling model so as to obtain wave elements, wave radiation stress and bed surface shear stress under the wave current coupling condition.
The model comprises a water flow module considering wave influence, a wave field module considering the water flow influence, and a transmission process of influence parameters of a wave field and a flow field.
Further, the step S2 includes the following steps:
s21, calculating the linear distribution of the wave breaking turbulence intensity in advance, and adopting an empirical linear vortex-viscosity coefficient formula in the invention:
ν t =2f w h( ) 1/3 ；
in the formula,f w for turbulent transformation coefficient, it can be 0.1;D f the energy lost by bottom friction is expressed asD f =1/(2π0.5)ρf w u 3 ord ，f w The coefficient of the wave friction resistance is the coefficient,u ord the maximum oscillation speed of the waves near the bottom;D r the energy lost by water roll is expressed asD r =2gEr sinβ/c，βIs the inclination angle of the water roll.
S22, the mass conservation equation after the waves are broken is established by the on-way distribution of the wave energy and the water roll energy provided by the S12 wave module:
in the formula,h t the terms 2 and 3 represent the time-average flow of the wave nonlinear motion and the wave surface water tumble respectively, which are the water depth at the wave trough.
S23, establishing a relation between the turbulent vortex viscosity coefficient and the crushing energy loss:
；
in the formula,the time-averaged shear stress; v istIs the turbulent vortex viscosity coefficient; u is the time-average flow rate; ρ is the fluid density; z is the elevation from the bed.
S24, calculating the time-averaged shear stress in the middle layer and the boundary layer under the assumption that the beach cross section is closed, wherein the vertical distribution formula is as follows:
，δ≤z≤h t
， 0≤z≤δ
wherein δ is the boundary layer thickness, and the expression δ =0.09 (a/k)s)0.82ksWhereinAThe amplitude of the wave near-bottom water particle is shown;k s the rough height of the bed surface;h t the elevation, m, of the trough surface from the bottom bed; s、 brespectively, the time-averaged trough surface shear stress and the time-averaged bed surface shear stress, m2/s2。
S25, under the condition that flow velocity is applied to the bed surface and no slippage exists, the average shear stress at the given time on the wave valley surface is obtained, and the calculation formula is as follows:
in the formula,in order to average the water surface elevation,His the wave height.
And S26, taking the obtained time-average flow velocity profile and the vertical distribution of the time-average shear stress as input parameters, and solving the random distribution data of the bed surface shear stress and the flow velocity according to a bottom boundary layer control equation.
In the boundary layer, the pressure gradient is distributed along the water depth and is in a hyperbolic cosine relationship, and the pressure gradient is equivalent to and unchanged on the boundary layer; the control equation is:
；
；
where w is the vertical flow rate and the subscript ∞ represents the physics at the boundary on the boundary layer, with the rest being as before.
And the step S3 is used for calculating sediment transport (suspended transport and bed transport) under the action of near shore power and beach change caused by sediment movement, and realizes accurate quantitative description of interaction of a hydrodynamic process, the sediment transport process and the beach landform evolution. Further, the method comprises the following specific steps:
s31, the transport mode of the near-shore silt mainly includes bed load and suspended load (suspension of viscous silt and part of non-viscous silt). The total sand transport rate is the sum of the sand transport rate of the bed load and the sand transport rate of the suspended load, and is as follows:q=q b +q s 。
calculating the bed load sand transport rate, wherein the bed load sand transport rate when the wave flow acts together is calculated according to the following formula:
；
in the formula,q b for single width volume bed mass sand transport rate m3/(m s)，s=ρ s /ρ 0 Is the relative density of the silt, d50The median diameter of silt and the empirical coefficient αn、b；θ cw,m ，θ cw The average and maximum Shields numbers of the coexisting wave streams,θ cw,m =（θ2 c+θ2 w, m+2θ c θ w,m cosψ）0.5，θ cw =(θ2 c+θ2 w+2θ c θ w cosψ)0.5，ψis the angle between the wave and the water flow.θ cr For the critical start-up of the Shields parameters,θ cr =0.30/(1+1.2D)+0.055[1-exp(-0.020D*)],D*=d50[g(s-1/v2)]1/3。
calculating the suspended load sand transporting rate, wherein the suspended load sand transporting rate formula when the wave flow acts together is as follows:
；
in the formula,U m is the average flow rate;Sthe sand-carrying capacity is the coexistence of wave flow;nis the porosity of the silt; omegasThe settling speed of the silt is determined;f w is the wave friction coefficient;D b is the fragmentation wave energy dissipation rate;β 1 ，β 2 andβ 3 are empirical coefficients.
S32, establishing a suspended load sediment transport module, wherein a two-dimensional sediment continuous equation along the water depth is adopted:
in the formula,S C the average sand content is the vertical line;、is the turbulent diffusion coefficient;is a silt flushing item;、the source amount and the source sand content of a horizontal unit area;for the erosion and deposition function, other symbolic meanings are defined with the hydrodynamic model.
And S33, establishing a transport sediment transport module. The invention directly adopts a calculation formula of the moving transport rate per unit area under the combined action of wave flow.
S34, calculating sediment settlement and resuspension, and settling the sediment when the shear stress at the bottom of the bed surface is smaller than the critical shear stress required by the sediment settlement; when bed surface bottom shear stress is greater than the critical shear stress that silt needs and the critical shear rate of resuspension, silt resuspension:
sediment sedimentation:
and (3) resuspending the silt:and is andu≥u *cs
in the formula:critical shear stress required for silt settling;critical shear stress required for silt movement;u *cscritical shear rate for resuspension of silt
S35, establishing a small-scale beach section evolution model, and respectively calculating seabed changes caused by sediment transport changes, wherein the seabed changes are divided into transverse motions to cause beach changes, namely the beach section model; the longitudinal transport of sediment, namely the evolution model of the beach caused by the transport of sediment along the bank, is expressed as the change of the position of the shoreline. Considering the influence of the gradient of the bottom surface of the seabed, the method is based on the conservation of the mass of the sediment, namely a sediment transport continuous equation:
in the formula,q x、 q y total sand transport rateqIn the x, y direction components;are empirical coefficients.
And S4, establishing a large-scale shoreline change model by using the long-period control function set, wherein the large-scale shoreline change model is used for scale expansion in the small-scale shoal section model, and quantitative prediction of large-scale shoreline change is realized. Further, it comprises the following steps:
and S41, establishing a large-scale shoreline change model by using the long-period control function set (comprising a scale extension function for extending the real-time calculated small-scale effect process to a longer time scale and a correction function in model operation).
Specifically, the scale extension function forms a long series required by calculation of an internal power process which comprises runoff, tide, wave and the like and can completely describe regional beach section evolution, and can be generalized into N power cycles, wherein each cycle time is T, the ith power cycle time period is integrated with a seabed deformation equation, each item is subjected to Taylor expansion, the influence of a second-order item is ignored, and the short-term change of the seabed change in a single calculation time step is not large, so that the structure converts the seabed short-term section into the seabed long-term evolution.
Specifically, a correction function, namely, a seabed updating function is controlled, and in an iteration period, when the total seabed change amount is larger than a set threshold value min (epsilon, 0.2h), the seabed is updated, and the hydrodynamic force is recalculated; secondly, simulating the frequency of representative extreme hydrological conditions by using a control function of the extreme hydrological conditions, and giving the recurrence period of storm surge by using Gumbul distribution; the recurrence period for large floods is given by the pearson type iii distribution.
Step S5 adopts Godunov finite volume method on the non-structural grid, solves the large scale bank line change model, obtains the large scale bank line change value, further specifically includes the following steps:
s51, synchronously solving the water flow-sediment transport flux by adopting an explicit format, and combining a MUSCL-Hancock prediction-correction format to enable a numerical format to have a conservation property;
the calculation of the diffusion flux of the substance uses the divergence theorem, and the formula is as follows:
in the formula,andrespectively based on the water depth reconstruction values of the left and right side units of the interface;andare respectively variable atx-andy-the original gradient of the direction.
Calculating the convection flux of the material by adopting a windward format:
in the formula,indicating that the water flows from the interface left cell to the right cell, and vice versa.The concentration of the substance in the left unit;the concentration of the substance is the right unit.
The MUSCL-Hancock prediction-correction format is a unit center value used for variable reconstruction in a correction step by taking a result of the prediction step, and can ensure the time second-order precision of the calculation format.
In the prediction step, the dry unit takes the concentration value at the current moment as the result of the prediction step; the wet unit then calculates the material concentration prediction step using a control equation based on the original variables according to the following equation:
in the formula,andare respectively asxAndythe limiting slope of the direction.The calculation time step is consistent with the calculation time step of the hydrodynamic model.
In the calibration step, the substance concentration is updated from time t to time t + Δ t. Based on a conservation form material transport control equation, the unit average conservation vector is updated by the following formula:
in the formula, subscriptiAndkrespectively numbering grids and unit edges;Ωis the unit area;lis the length of the edge;is the convective flux;is the diffusion flux;S i is an approximation of the source term. The model adopts an explicit format to calculate a source and sink term and an attenuation term. When updating the concentration of the substance, if the unit water depth is less than a certain threshold (generally 10 is taken)-6m), the concentration is set to 0.
Acquiring a small-scale beach deformation numerical value, dispersing and solving a small-scale beach deformation equation,
in the formula,is a moving mass per unit area (kg/m)2)；iNumbering the units;knumbering the edges of the cells;sand transport flux for marginal bed load:
adopting a display format to process a riverbed total deformation equation:
in the formula,Z b is the unit central bottom elevation;grouping the number of the silt particle sizes;nthe silt particle size is numbered in groups.
Because the model bottom elevation is defined in the node, the node elevation needs to be updated in an interpolation mode after the model is updated and the unit center bottom elevation is updated.
Acquiring a large-scale shoreline deformation value, and performing dispersion and solution on a large-scale seabed deformation equation according to the following formula:
and the equations are integrated over each unit, the final discrete format can be written as:
in this patent, the silt transportation volume of small-scale beach section omits the very little factor of influence through the main dynamic effect frequency of generalization, constitutes the net influence of long-time sequence, expands the transportation volume of large-scale bank line and corresponding bank line change, realizes long-time scale power geomorphology change under the limited power calculation time from this.
As the landform is updated once after each landform time step, next step of hydrodynamic calculation is based on the updated landform, and the operation is repeated, the interaction and feedback between the hydrodynamic force and the landform can be well considered, and the situation accords with the actual landform evolution process situation.
Compared with the prior art, the invention has the following advantages and effects:
the invention provides a numerical simulation method for coupling small-scale beach section and large-scale shoreline change based on a sand conveying process, which realizes accurate quantitative description of interaction of a hydrodynamic process, a sediment transport process and landform evolution, omits factors with little influence through generalization of main dynamic action frequency, forms net influence of a long-time sequence, realizes expansion of transport capacity from sediment transport capacity on a small-scale beach to transport capacity on a large-scale shoreline landform, solves the problems of matching of multiple time-space scales in evolution, calculation time consumption in long-period simulation and the like, and remarkably improves the practical application value of beach evolution simulation.
Drawings
FIG. 1 is a flow chart of the method of the present invention.
Fig. 2 is a schematic diagram of a model range of the section of the kahn bay according to the embodiment of the present invention.
Fig. 3 is a graph of the effective wave height of the ordinary wave direction-E toward wave purification bay.
Fig. 4 is a schematic view of the normal wave direction-E to the wave purification bay wave generation field.
Fig. 5 is a schematic diagram of the direction of the sand transportation along the bank in the Jing gulf.
Fig. 6 is a diagram of the annual scour rate profile for the Jing gulf obtained by the present invention.
FIG. 7 is a graph showing the contour line changes at different times.
Fig. 8 is a shoreline verification diagram.
Fig. 9 is a comparison graph of the position of the shore line in the state of equilibrium in the Jing bay and the position of the shore line fitted by an empirical formula.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention clearer and clearer, the present invention is further described in detail through embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.
FIG. 1 is a flow chart of the method of the present invention. This patent generally includes the steps of:
s1, establishing a wave flow coupling model for accurately and quantitatively describing the plane near-shore power;
s2, establishing a bottom offshore flow model for accurately and quantitatively describing the offshore bottom offshore flow;
s3, establishing a small-scale beach profile evolution model;
s4, establishing a large-scale shoreline change model by using the long-period control function set;
and S5, solving the large-scale shoreline change model by adopting a Godunov finite volume method on the unstructured grid, and acquiring a large-scale shoreline change value.
The embodiment of the invention is mainly explained by carrying out simulation demonstration on the sediment transport condition of a Jing gulf-sandy coast.
Jing bay is located the southeast seaside in Huilai county, and the last cape angle is north fort platform, and the lower cape angle is the position angle, and the gulf mouth opens to the east, is the typical sand matter cape room arc bay system in the east coast of Guangdong, and the geographical position sketch map is shown in FIG. 1. After the breakwater project of a certain power plant is finished, the power and deposition process in the bay are changed, and the shape of the beach is also remodeled. The 'sink' and 'source' of the coastal sand transportation are reversed, and the relative erosion or erosion backing-up of the oasis bank section (the straight tangent line of the original arc-shaped coast system) in the bay is accelerated.
The movement form of the tidal current in the Jing gulf sea area is generally a reciprocating flow in the northeast-southwest direction, and the directions of the rising current and the falling current are consistent with the trends of the isobathes. The current flow velocity outside the extragulal-20 m contour is relatively large, while the current flow effect inside the gulf is weak. Therefore, the influence of the tidal current on the movement of the near-shore silt is relatively small.
The yearly distribution statistics of each wave level wave at the depth of 20m in Jing Bay are shown in Table 1. As can be seen from table 1, the waves in the jing gulf are in the direction of frequent waves, i.e., the dominant waves are in the direction of E, and the annual frequency of occurrence is 18.9%; the minor wave direction is NE direction and SE direction, and the annual frequency of occurrence is 17.8% and 17.2%, respectively. The frequency of wave occurrences between the NE direction, centered at the E direction, and the SSE direction, accounts for 82.5% of the annual occurrence rate, with most waves occurring in the E direction-NE direction and most swells occurring in the SE direction-SSE direction.
The content of average suspended sand in the bay of Jing Hai in summer is 0.0117kg/m3Maximum value of 0.1384kg/m3. Under the influence of tidal current and wave action, the distribution of the suspended sand content in the bay has the distribution characteristic that the suspended sand content is increased and then reduced from the bank to the open sea. In winter, the average sand content in bay is 0.0286kg/m due to the disturbance of strong northeast monsoon3Maximum value of 0.1501kg/m3There is a tendency to increase from the top layer to the bottom layer.
Most of the sediments from the beach of Jing gulf to the middle of the gulf are medium sands. The sand in the beach has the largest median particle size and becomes finer towards the deep sea. Wherein the median diameter of the sediment at the shoal of the coast of the capital section is 0.27mm, the average median diameter of the sediment within the depth line of minus 5m is 0.22mm, and the average median diameter of the sediment between the depth lines of minus 10m to minus 15m is 0.14mm to 0.17 mm.
TABLE 1 annual wave distribution at the depth of 20m in Jing Bay
| Direction of rotation | Frequency (%) | Period of time | H10% |
| N | 0.20 | 5.70 | 2.20 |
| NNE | 6.50 | 5.3 | 1.91 |
| NE | 17.80 | 5.30 | 2.07 |
| ENE | 8.40 | 5.00 | 1.82 |
| E | 18.90 | 6.50 | 1.49 |
| ESE | 3.80 | 6.40 | 1.00 |
| SE | 17.20 | 6.40 | 1.08 |
| SSE | 16.40 | 5.90 | 1.28 |
| S | 9.10 | 6.40 | 1.00 |
| SSW | 0.50 | 5.70 | 1.66 |
| SW | 0.80 | 6.80 | 0.95 |
| WSW | 0.00 | 0.00 | 0.00 |
| W | 0.00 | 0.00 | 0.00 |
| WNW | 0.00 | 0.00 | 0.00 |
| NW | 0.10 | 5.40 | 2.20 |
| NNW | 0.10 | 5.40 | 2.20 |
The model range in the embodiment of the invention is shown in FIG. 2, and comprises an arc coast system between the cape corners of Jinghai bay to a deep line of-20 m outside the bay, the length of the model is 7.8km in south and north, the width of the model is 4.8km in east and west, and the area of the simulated water area is about 35km2The model also comprises a section from the estuary of the lagoon to the water gate of the lion-stone lake.
The computational mesh is a triangular mesh, and has 84074 cells, 28394 nodes and 55680 edges in total. The minimum length of the grid edge is 2m, and the minimum area of the grid is 5m2。
And the north boundary of the model boundary is positioned at the north part of the north fort and is about 22km away from a Haimen station at the north part, and the harmonic constant of the tide level survey station is used for harmonic analysis to obtain a tide level process in the calculation time period. And (4) harmonic analysis is carried out on the south boundary by utilizing a harmonic constant of the station survey station of the harbor A to obtain a tide level course in the calculation time period. The tide levels for the east boundary and other water points are a linear interpolation of the tide levels for the two points.
In order to realize the long-term evolution of the seabed and avoid overlong calculation time, the invention adopts a reduction technology of power and terrain coupling to generalize a main power process influencing the evolution of the seabed.
When a long-time sequence is constructed, a main power process influencing the evolution of the long seabed needs to be identified and used as a main power factor to replace an actual sequence power process. The average runoff and the sand transportation amount of years are respectively used for replacing the sand transportation amount of rivers, and in order to consider the influence of flood, design flood with different representative frequencies is added; the generalized mode of 'representing tide' is used for representing the long-period average transport rate which can reflect the actual tide in the long-term tide data; the general wave height capable of reflecting short-period extreme wave height and long-period trend beach erosion in the long-term wave height statistical data is represented by a 'representative wave' or 'characteristic wave' generalized mode, weighted average generalized processing is carried out according to the frequencies of two different wave conditions, and meanwhile, the influence of the occurrence frequency of storm surge is considered.
As can be seen from the analysis of the power of the near shore region of the jing gulf, the main driving force for the transport of silt in the jing gulf is the wave and wave-borne coastal currents. Considering the shoreline has the NE-SW trend, the wave directions that play a major role in transporting sand along the shore in the bay are the ENE, E, ESE, SE, SSE and S directions. Effective wave height given on open sea-20 m equal depth lineH s =H10%/1.27，H10%The square root of each wave stage is taken and the simulated time of each wave is converted from the respective frequency, as detailed in table 2.
Sand content value given on open borders: when the sand-carrying amount of the boundary is given, considering that the coast erosion-silting amount caused by the tidal current scouring is small under the normal condition, in order to observe the more remarkable change condition of the coasts between the promontory, the sand content on the boundary is equal to the sand-carrying force of the water body when the boundary is arranged in inflow, namely the sand content of the water body under the saturated state.
The median grain diameter of bed sand is set according to the isobathic line differentiation, and the median grain diameter of silt in the low-lying cape angle shoal scouring area of nearly bank is 0.27mm promptly, and the median grain diameter of silt in the upper cape angle siltation area is 0.22mm, and the median grain diameter of silt gradually passes through 0.22mm to 0.17mm between the isobathic lines of 5m ~10m, and the median grain diameter of silt gradually passes through 0.17mm to 0.14mm between the isobathic lines of 10m ~20 m.
TABLE 2 typical wave height (m), period(s) and time length (day) of each level of waves in the structural wave series
| H s | ENE | E | ESE | SE | SSE | S |
| 0.22m | 0 | 0 | 0.4 | 1.1 | 0.7 | 0.7 |
| 0.83m | 10.6 | 35.8 | 13.1 | 60.6 | 43.4 | 30.7 |
| 1.82m | 19.0 | 27.0 | 0.7 | 7.7 | 15.3 | 1.8 |
| 3.20m | 1.1 | 0 | 0 | 0 | 0.4 | 0 |
| Period(s) | 5 | 6.5 | 6.4 | 6.4 | 5.9 | 6.4 |
In this embodiment, the model predicts the shoreline in 2015 as verification and predicts the future shoreline variation with the actual terrain measured in the top bay of 2008 as the initial terrain and initial shoreline position. The method comprises the following specific steps:
step 1, a wave-current coupling model is established firstly, and is used for calculating main hydrodynamic factors (such as tidal current, wind-driven circulation, runoff, near-shore current caused by wave breaking and the like) influencing the beach evolution, so that accurate quantitative description of plane near-shore power is realized.
Firstly, establishing a two-dimensional hydrodynamic module: in the invention, a conservation-form two-dimensional shallow water equation is used as a water flow control equation, which is as follows:
in the formula,is a conservation vector;、are respectively as、A directional convective flux vector;、are respectively as、Diffusion flux vector caused by directional reynolds stress;、are respectively as、A diffusion flux vector caused by directional secondary flow;for the source term vector:
in the formula,is the water depth;、respectively, average flow velocity in the vertical direction、A component of direction;is the bottom elevation;the intensity of rainfall is used;the infiltration strength is obtained;is the turbulent viscosity coefficient in the horizontal direction,，is a coefficient of proportionality that is,in order to be the karman coefficient,the flow rate is cut for the bed surface;、、、is a diffusion stress term caused by secondary flow;is the acceleration of gravity;is the coefficient of the Coriolis force,，is the angular velocity of the earth's rotation,the local latitude is;in order to be the wind stress,、，the water surface wind stress drag coefficient;、density of air and water, respectively;、are respectively as、A wind speed component (m/s) at a height of 10m above the water surface in the direction;、、the wave radiation stress is provided by a wave model;、are respectively as、The slope of the friction resistance in the direction,、. Wherein,h、u、vis the amount of strain.
Then, a wave module is established, and the invention adopts the third generation wave of SWAN based on the conservation equation of wave action quantity. A wave field module for providing wave elementsE w 、E r 、、HAnd calculating wave radiation stress and wave surface water rolling momentum flow. The specific formula is as follows:
in the formula:is the wave propagation direction andthe included angle of the axes;is the average wave energy per wave period of the water column,，nis the ratio of the group velocity to the phase velocity,，in terms of the wave number, the number of waves,E r is the energy density of the water roll.
Because the wave radiation stress changes along the way, the mean water surface changes, and the mean water surface elevation (water increase and decrease) can be calculated by the following formula in combination with the consideration of the influence of the rolling water:
the model comprises a water flow module considering wave influence, a wave field module considering the water flow influence, and a transmission process of influence parameters of a wave field and a flow field.
The results of analyzing the wave-generated flow field are represented by the direction from the open wave to the direction-E in Jing gulf. The map of E in Jing Bay to the height of the significant wave is shown in FIG. 3. After the dominant wave enters the bay, the dominant wave is deformed under the influence of the terrain and the breakwater of the power plant, so that the wave force in the bay is redistributed and gradually deflected in the bay, and finally the situation that the dominant wave intersects with the trend of the shoreline at an acute angle is achieved, and the influence of the dominant wave on the northern section of the shoreline is greatly reduced. At this time, the average wave height in the south of the bay is about 0.9 to 1.2m, and the average wave height in the north of the bay is very small under the protection of the breakwater of the power plant.
The results of the modeling of the normal wave direction-E to wave-borne coastal currents in the Jing Bay are shown in fig. 4. The gulf top region has a relatively low flow velocity, and the south region of the gulf is mainly coastal flow from south to north. There are more complex circulation systems in the southern, nearshore region, but the flow velocities are also relatively small, which should be associated with a relatively irregular wave height distribution in this region.
Step 2, establishing a bottom offshore flow model; the wave element provided in the step 1 is used as an input condition, the turbulence intensity is expressed based on the linear vortex viscosity coefficient vertical distribution, the bottom offshore flow is driven by the wave residual momentum flow gradient and the pressure gradient caused by increasing and decreasing the water level together, a bottom offshore flow model is established, and the bottom offshore flow model is used for calculating the bed shear stress and the water flow velocity distribution at any time under the coupling action of waves and flow, so that the accurate quantitative tracing of the near-shore bottom offshore flow is realized.
Firstly, linear distribution expressing wave-breaking turbulence intensity is calculated in advance:
；
in the formula,f w for turbulent transformation coefficient, it can be 0.1;D f the energy lost by bottom friction is expressed asD f =1/(2π0.5)ρf w u 3 ord ，f w The coefficient of the wave friction resistance is the coefficient,u ord the maximum oscillation speed of the waves near the bottom;D r the energy lost by water roll is expressed asD r =2gEr sinβ/c，βIs the inclination angle of the water roll.
Secondly, establishing a mass conservation equation after the waves are broken:
in the formula,h t the terms 2 and 3 represent the time-average flow of the wave nonlinear motion and the wave surface water tumble respectively, which are the water depth at the wave trough.
Thirdly, establishing a relation between the turbulent vortex viscosity coefficient and the crushing energy loss:
；
in the formula,the time-averaged shear stress; v istIs the turbulent vortex viscosity coefficient; u is the time-average flow rate; ρ is the fluid density; z is the elevation from the bed.
Under the assumption that the beach cross section is closed, the time-averaged shear stress in the middle layer and the boundary layer is calculated, and the vertical distribution formula is as follows:
，δ≤z≤h t
， 0≤z≤δ
wherein δ is the boundary layer thickness, and the expression δ =0.09 (a/k)s)0.82ksWhereinAThe amplitude of the wave near-bottom water particle is shown;k s the rough height of the bed surface;h t the elevation, m, of the trough surface from the bottom bed; s、 brespectively, the time-averaged trough surface shear stress and the time-averaged bed surface shear stress, m2/s2。
Under the condition that the flow velocity is applied to the bed surface and no slippage exists, the time-averaged shear stress on the wave valley surface is obtained, and the calculation formula is as follows:
；
in the formula,in order to average the water surface elevation,His the wave height.
And finally, taking the obtained time-average flow velocity profile and the vertical distribution of the time-average shear stress as input parameters, and solving the random distribution data of the bed surface shear stress and the flow velocity according to a bottom boundary layer control equation.
The control equation is:
where w is the vertical flow rate and the subscript ∞ represents the physics at the boundary on the boundary layer, with the rest being as before.
And 3, establishing a small-scale beach profile evolution model for calculating sediment transport (suspended transport and bed transport) under the action of near shore power and beach change caused by sediment movement, and realizing accurate quantitative description of interaction of a hydrodynamic process, a sediment transport process and beach landform evolution.
The transport pattern of near-shore silt mainly includes bed load and suspended load (suspension of viscous silt and part of non-viscous silt). The total sand transport rate is the sum of the sand transport rate of the bed load and the sand transport rate of the suspended load, and is as follows:q=q b +q s 。
calculating the bed load sand transport rate, wherein the bed load sand transport rate when the wave flow acts together is calculated according to the following formula:
in the formula,q b for single width volume bed mass sand transport rate m3/(ms)，s=ρ s /ρ 0 Is the relative density of the silt, d50The median diameter of silt and the empirical coefficient αn、b；θ cw,m ，θ cw The average and maximum Shields numbers of the coexisting wave streams,θ cw,m =（θ2 c+θ2 w, m+2θ c θ w,m cosψ）0.5，θ cw =(θ2 c+θ2 w+2θ c θ w cosψ)0.5，ψis the angle between the wave and the water flow.θ cr For the critical start-up of the Shields parameters,θ cr =0.30/(1+1.2D)+0.055[1-exp(-0.020D*)],D*=d50[g(s-1/v2)]1/3。
calculating the suspended load sand transporting rate, wherein the suspended load sand transporting rate formula when the wave flow acts together is as follows:
in the formula,U m is the average flow rate;Sthe sand-carrying capacity is the coexistence of wave flow;nis the porosity of the silt; omegasThe settling speed of the silt is determined;f w is the wave friction coefficient;D b is the fragmentation wave energy dissipation rate;β 1 ，β 2 andβ 3 are empirical coefficients.
In the invention, a suspended sediment transport module is established by adopting a two-dimensional sediment continuous equation along the water depth:
in the formula,S C the average sand content is the vertical line;、is the turbulent diffusion coefficient;is a silt flushing item;、the source amount and the source sand content of a horizontal unit area;for the erosion and deposition function, other symbolic meanings are defined with the hydrodynamic model.
In the invention, a calculation formula of the moving transport rate per unit area under the combined action of wave flow is directly adopted to establish a moving sediment transport module.
And (3) calculating the sand transportation along the coast of the Jing Bay by using the built model, wherein the calculation result is shown in a table 3, and a schematic diagram of the section arrangement and the sand transportation direction is shown in a figure 5. The coastal sediment transporting direction of the northern arc bank section (I-II section) of the Jing bay points to the northern gun platform, namely the gulf top. The annual sand transportation amount is 5.36-17.52 ten thousand meters3On the left and right, the situation is consistent with that the shoreline of 2008-2015 is pushed to the outer sea, i.e. the siltation and the forward spreading. The sand conveying direction along the bank of the arc coast straight-line segment (III-IV section) of the Jing Bay points to the depth angle, and the sand conveying amount is 2.18-3.56 ten-thousand meters3On the left and right sides, because the direction of sand conveying along the bank is opposite to the direction of sand conveying along the bank on the arc top side, the supply of the sand flow along the bank cannot be obtained, and the erosion of the bank line is shown. The coastal sand conveying direction of the arc south bank segment (V-VIII section) of the Jing Bay points to the northern gun platform, and the sand conveying amount is 17.89-52.11 ten-thousand meters3Left and right. Based on the consideration that the relative independence of the Jing gulf and the long-term shortage of sand coming from the sea area (inland frame), the erosion retreat of the coast section will be a long-term continuous process and trend.
TABLE 3 results of coastal sediment transport
Then, sediment settlement and resuspension are calculated, and when the shear stress of the bottom of the bed surface is smaller than the critical shear stress required by sediment settlement, the sediment settlement is carried out; when bed surface bottom shear stress is greater than the critical shear stress that silt needs and the critical shear rate of resuspension, silt resuspension:
sediment sedimentation:；
and (3) resuspending the silt:and is andu≥u *cs(ii) a In the formula:critical shear stress required for silt settling;critical shear stress required for silt movement;u *csthe critical shear rate for resuspension of silt.
Establishing a small-scale beach section evolution model, respectively calculating seabed changes caused by sediment transport changes, and dividing the seabed changes into transverse motions to cause beach changes, namely the beach section model; the longitudinal transport of sediment, namely the evolution model of the beach caused by the transport of sediment along the bank, is expressed as the change of the position of the shoreline. Considering the influence of the gradient of the bottom surface of the seabed, the method is based on the conservation of the mass of the sediment, namely a sediment transport continuous equation:
in the formula,q x、 q y total sand transport rateqIn the x, y direction components;are empirical coefficients.
And 4, expanding the scale in the small-scale beach section model by using a long-period control function set (comprising a scale extension function for extending the real-time calculated small-scale effect process to a longer time scale and a correction function in model operation) so as to realize quantitative prediction of large-scale shoreline change.
Specifically, the scale extension function forms a long series required by calculation of an internal power process which comprises runoff, tide, wave and the like and can completely describe regional beach section evolution, and can be generalized into N power cycles, wherein each cycle time is T, the ith power cycle time period is integrated with a seabed deformation equation, each item is subjected to Taylor expansion, the influence of a second-order item is ignored, and the short-term change of the seabed change in a single calculation time step is not large, so that the structure converts the seabed short-term section into the seabed long-term evolution.
Specifically, a correction function, namely, a seabed updating function is controlled, and in an iteration period, when the total seabed change amount is larger than a set threshold value min (epsilon, 0.2h), the seabed is updated, and the hydrodynamic force is recalculated; secondly, simulating the frequency of representative extreme hydrological conditions by using a control function of the extreme hydrological conditions, and giving the recurrence period of storm surge by using Gumbul distribution; the recurrence period for large floods is given by the pearson type iii distribution.
And 5, solving a large-scale shoreline change model by adopting a Godunov finite volume method on the unstructured grid to obtain a large-scale shoreline change numerical value.
The method specifically comprises the following steps: synchronously solving the water flow-sediment transport flux by using an explicit format, and combining MUSCL-Hancock prediction-correction;
the calculation of the diffusion flux of the substance uses the divergence theorem, and the formula is as follows:
in the formula,andrespectively based on the water depth reconstruction values of the left and right side units of the interface;andare respectively variable atx-andy-the original gradient of the direction.
Calculating the convection flux of the material by adopting a windward format:
in the formula,indicating that the water flows from the interface left cell to the right cell, and vice versa.The concentration of the substance in the left unit;the concentration of the substance is the right unit.
In the prediction step, the dry unit takes the concentration value at the current moment as the result of the prediction step; the wet unit then calculates the material concentration prediction step using a control equation based on the original variables according to the following equation:
in the formula,andare respectively asxAndythe limiting slope of the direction.The calculation time step is consistent with the calculation time step of the hydrodynamic model.
In the calibration step, the substance concentration is updated from time t to time t + Δ t. Based on a conservation form material transport control equation, the unit average conservation vector is updated by the following formula:
in the formula, subscriptiAndkrespectively numbering grids and unit edges;Ωis the unit area;lis the length of the edge;is the convective flux;is the diffusion flux;S i is an approximation of the source term. The model adopts an explicit format to calculate a source and sink term and an attenuation term. When updating the concentration of the substance, if the unit water depth is less than a certain threshold (generally 10 is taken)-6m), the concentration is set to 0.
Obtaining the small-scale beach deformation value, dispersing and solving the small-scale beach deformation equation,
in the formula,is a moving mass per unit area (kg/m)2)；iNumbering the units;knumbering the edges of the cells;sand transport flux for marginal bed load:
adopting a display format to process a riverbed total deformation equation:
in the formula,Z b is the unit central bottom elevation;grouping the number of the silt particle sizes;nthe silt particle size is numbered in groups.
Because the model bottom elevation is defined in the node, the node elevation needs to be updated in an interpolation mode after the model is updated and the unit center bottom elevation is updated.
And (3) simulating 1 year by using a constructed 'representative wave' of 1 year, and calculating the scouring and silting conditions of the Jing gulf quayside, wherein the scouring and silting distribution is shown in a figure 6. As can be seen from fig. 6, the beach is eroded at the promontory corner under the bay and silted at the crest of the bay, i.e. the breakwater shelter, which is consistent with the conclusion of the terrain evolution (fig. 7 actually measures the contour change map at different times). The sedimentation rate of the wave-proof bank shielding section is more than 0.2m/a, and the scouring rate of the lower cape corner section is more than 0.1 m/a.
Obtaining a large-scale shoreline deformation value, and carrying out dispersion and solution on a large-scale seabed deformation equation according to the following formula:
；
and the equations are integrated over each unit, the final discrete format can be written as:
。
firstly, expanding the constructed 'representative wave' of 1 year to more than 10 years by using a long-period control function set; the large-scale shoreline deformation model is constructed by the method, and the model evolves to the shoreline of 2015 on the basis of the shoreline of 2008. The predicted 2015-year shore line was verified against the actually measured 2015-year shore line. The verification results are shown in FIG. 8. It can be seen that the calculated maximum flush point for 2015 substantially coincides with the actual measurement. Calculations show that the selection of relevant parameters in the model is reasonable, and the model can be used to predict the scour changes of the Jing Bay quayside.
On the basis of actually measuring the terrain and the shoreline in 2015, constructing a large-scale shoreline deformation model to predict the change condition of the future shoreline by using the method; when the net sand conveying amount on the near shore is not changed greatly, the balance profile is considered to be reached, the shore line is extracted, and the shore line is matched with a parabolic empirical formula to obtain the shore line for comparison, which is shown in figure 9. The difference between the model forecast equilibrium state and the fitting result of the empirical formula is small. Therefore, the model can reasonably reflect the erosion and deposition of the coast of Jing gulf and accord with the long-term receding and short-term evolution rule of the coastline.
In summary, the present embodiment is summarized as follows:
the method aims at the problems of multiple space-time scales, complex coupling effect of different modules (power, sediment transport and terrain evolution), expansion of numerical error accumulation effect on a long-time scale and the like existing in beach evolution. Establishing a wave flow coupling model, a bottom offshore flow model, a sediment transport and terrain change model, and realizing accurate and quantitative description of a physical mechanism of small-scale beach profile change caused by bottom offshore flow, wave flow interaction, additional shear stress of a bottom boundary layer, vertical and horizontal sediment transport and the like; the method is characterized in that main power is generalized, a long-period control function set is introduced, the sediment transport process is corrected based on water quantity and sediment conservation, a beach terrain change model based on longitudinal and transverse sediment transport rates is established, the expansion of the transport volume from the sediment transport volume of hydrodynamic scales to the transport volume of landform scales is realized, the simulation of large-scale shoreline change is completed, the problem of matching of multiple space-time scales in evolution is solved, and the fact that the calculation time in long-period simulation is not long is guaranteed. The method can obviously improve the practical application value of the beach evolution simulation.
The above examples are preferred embodiments of the present invention, but the present invention is not limited to the above examples, and any other changes, modifications, substitutions, combinations, and simplifications which do not depart from the spirit and principle of the present invention should be construed as equivalents thereof, and all such changes, modifications, substitutions, combinations, and simplifications are intended to be included in the scope of the present invention.
Claims (6)
1. A coupling simulation method for small-scale beach profile and large-scale shoreline change is characterized by comprising the following steps:
s1, establishing a wave flow coupling model for accurately and quantitatively describing the plane near-shore power;
s2, establishing a bottom offshore flow model for accurately and quantitatively describing the offshore bottom offshore flow;
s3, establishing a small-scale beach profile evolution model;
s4, establishing a large-scale shoreline change model by using the long-period control function set;
and S5, solving the large-scale shoreline change model by adopting a Godunov finite volume method on the unstructured grid, and acquiring a large-scale shoreline change value.
2. The coupling simulation method according to claim 1, wherein the step S1 comprises the steps of:
s11, establishing a hydrodynamic module, wherein the equation is as follows:
in the formula,is a conservation vector;、are respectively as、A directional convective flux vector;、are respectively as、Caused by directional Reynolds stressesA diffusion flux vector;、are respectively as、A diffusion flux vector caused by directional secondary flow;is a source term vector;
；；；
；；
；；
；
in the formula,is the water depth;、respectively, average flow velocity in the vertical direction、A component of direction;is the bottom elevation;the intensity of rainfall is used;the infiltration strength is obtained;is the turbulent viscosity coefficient in the horizontal direction,，is a coefficient of proportionality that is,in order to be the karman coefficient,the flow rate is cut for the bed surface;、、、is a diffusion stress term induced by the secondary flow;is the acceleration of gravity;is the coefficient of the Coriolis force,，is the angular velocity of the earth's rotation,the local latitude is;in order to be the wind stress,、，the water surface wind stress drag coefficient;、density of air and water, respectively;、are respectively as、A wind speed component at a height of 10m above the water surface in the direction;、、the wave radiation stress is provided by a wave model;、are respectively as、The slope of the friction resistance in the direction,、(ii) a Wherein,h、u、vis a dependent variable;
s12, establishing a wave module, wherein the equation is as follows:
；
；
in the formula:is the wave propagation direction andthe included angle of the axes;is the average wave energy per wave period of the water column,，nis the ratio of the group velocity to the phase velocity,，in terms of the wave number, the number of waves,E r is the energy density of the water roll;
and S13, coupling the hydrodynamic module in S11 with the wave module in S12 to form a wave current coupling model.
3. The coupling simulation method according to claim 1, wherein the step S2 comprises the steps of:
s21, setting the linear distribution of wave-breaking turbulence intensity:
；
in the formula,f w the turbulent transformation coefficient;D f the energy lost for the friction of the bottom part,D f =1/(2π0.5)ρf w u 3 ord ，f w the coefficient of the wave friction resistance is the coefficient,u ord the maximum oscillation speed of the waves near the bottom;D r in order to obtain the energy lost by the water roll,D r =2gEr sinβ/c，βthe inclination angle of the water roll;
s22, setting a mass conservation equation after the waves are broken:
in the formula,h t the items 2 and 3 represent the time-average flow of wave nonlinear motion and wave surface water tumble respectively, wherein the water depth at the wave trough is the water depth;
s23, setting a correlation formula of the turbulent vortex viscosity coefficient and the crushing energy loss:
；
in the formula,the time-averaged shear stress;ν t is the turbulent vortex viscosity coefficient;is the time-average flow rate;ρis the fluid density; z is the elevation from the bed;
s24, calculating the time-averaged shear stress in the middle layer and the boundary layer under the condition that the beach cross section is set as a closed condition, wherein the vertical distribution formula is as follows:
，δ≤z≤h t
，0≤z≤δ
in the formula,δin order to be the boundary layer thickness,δ=0.09( A /ks)0.82kswhereinAThe amplitude of the wave near-bottom water particle is shown;k s the rough height of the bed surface;h t the elevation of the trough surface from the bottom bed; s、 brespectively time-averaged trough surface shear stress and time-averaged bed surface shear stress;
s25, under the condition that no slippage is applied to the set bed surface, the average time shear stress on the wave valley surface is obtained, and the calculation formula is as follows:
；
in the formula,in order to average the water surface elevation,His the wave height; the time-average water surface elevation calculation formula is as follows:
；
s26, coupling the parameters in the previous steps, and solving the data of the shearing stress and the flow velocity of the bed surface distributed at any time, wherein the equation is as follows:
；
；
in the formula,wfor vertical flow rates, the subscript ∞ represents the physics at the upper boundary of the boundary layer.
4. The coupling simulation method according to claim 1, wherein the step S3 comprises the steps of:
s31, calculating the total sand transporting rate when the wave flows act togetherq，q=q b +q s ：
Wherein the bed load sand transport rateq b The calculation formula of (a) is as follows:；
in the formula,q b for single width volume bed mass transport rate, s =ρ s /ρ 0 Is the relative density of the silt,d 50 the median diameter of the sediment;θ cw,m ，θ cw the average and maximum Shields numbers of the coexisting wave streams,θ cw,m =（θ2 c+θ2 w,m+2θ c θ w,m cosψ）0.5，θ cw =(θ2 c+θ2 w+2θ c θ w cosψ)0.5，ψis the included angle between the wave and the water flow;θ cr for the critical start-up of the Shields parameters,θ cr =0.30/(1+1.2D)+0.055[1-exp(-0.020D*)]，D*=d50[g(s-1/v2)]1/3；
suspended load sand transport rateq s The calculation formula of (a) is as follows:
；
in the formula,U m is the average flow rate;Sc * the sand-carrying capacity is the coexistence of wave flow;nis the porosity of the silt; omegasThe settling speed of the silt is determined;f w is the wave friction coefficient;D b is the fragmentation wave energy dissipation rate;
s32, establishing a suspended load sediment transport module, wherein the equation is as follows:
in the formula,S C the average sand content is the vertical line;、is the turbulent diffusion coefficient;is a silt flushing item;、the source amount and the source sand content of a horizontal unit area;is a erosion and deposition function;
s33, calculating sediment settlement and resuspension:
sediment sedimentation:；
and (3) resuspending the silt:and is andu≥u *cs；
in the formula:critical shear stress required for silt settling;critical shear stress required for silt movement;u *csthe critical shearing speed for the resuspension of the silt;
s34, establishing a small-scale beach section evolution model, wherein the equation is as follows:
in the formula,q x、 q y total sand transport rateqIn the x, y direction components;are empirical coefficients.
5. The coupling simulation method of claim 1, wherein the set of long-period control functions in step S4 includes a scale extension function and a correction function; the correction function comprises a control seabed update function and an extreme hydrological condition control function: the scale extension function formula is as follows:
。
6. the coupling simulation method according to claim 1, wherein the step S5 comprises the steps of:
s51, synchronously solving the water flow-sediment transport flux by adopting an explicit format, and combining MUSCL-Hancock prediction-correction;
wherein, the calculation formula of the diffusion flux of the substance is as follows:
in the formula,andrespectively based on the water depth reconstruction values of the left and right side units of the interface;andare respectively variable atx-andy-an original gradient of direction;
the formula for calculating the convection flux of a substance is as follows: :
in the formula,indicating that the water flow is from the interface left cell to the right cell, whereas indicating that the water flow is from the interface right cell to the left cell,
the concentration of the substance in the left unit;the concentration of the substance in the right unit;
in the prediction step, the dry unit takes the concentration value at the current moment as the result of the prediction step; the wet unit calculates the concentration of the substance by adopting a control equation based on an original variable; the formula is as follows:
in the formula,andare respectively asxAndya limiting slope of direction;the calculation time step is consistent with the calculation time step of the hydrodynamic model;
in the correction step, the concentration of the substance is updated from the time t to the time t + delta t; based on a conservation form material transport control equation, updating the cell average conservation vector by the following formula:
in the formula, subscriptiAndkare respectively a gridNumbering and unit edge numbering;Ωis the unit area;lis the length of the edge;is the convective flux;is the diffusion flux;S i is an approximation of the source term; the model adopts an explicit format to calculate a source and sink item and an attenuation item; when updating the concentration of the substance, if the unit water depth is less than a certain threshold (generally 10 is taken)-6m), the concentration is set to 0;
s52, obtaining the small-scale beach deformation value, dispersing and solving the small-scale beach deformation equation,
in the formula,is the advancing mass per unit area;inumbering the units;knumbering the edges of the cells;for the marginal bed load sand transport flux, the formula is as follows:
adopting a display format to process a riverbed total deformation equation, wherein the formula is as follows:
in the formula,Z b is the unit central bottom elevation;grouping the number of the silt particle sizes;nnumbering the silt particle size groups;
s53, obtaining a large-scale shoreline deformation value, and solving the discretization and solving formula of the large-scale seabed deformation equation as follows:
and integrating the equation on each unit, the final discrete format formula is as follows:
。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810944700.4A CN109582996A (en) | 2018-08-19 | 2018-08-19 | A kind of coupled simulation method of small scale beach profile and large scale Shoreline changes |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810944700.4A CN109582996A (en) | 2018-08-19 | 2018-08-19 | A kind of coupled simulation method of small scale beach profile and large scale Shoreline changes |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN109582996A true CN109582996A (en) | 2019-04-05 |
Family
ID=65919660
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810944700.4A Pending CN109582996A (en) | 2018-08-19 | 2018-08-19 | A kind of coupled simulation method of small scale beach profile and large scale Shoreline changes |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109582996A (en) |
Cited By (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110388899A (en) * | 2019-06-04 | 2019-10-29 | 浙江省水利河口研究院 | The personal module calculated for tidal bore flow velocity vertical characteristics |
| CN110412312A (en) * | 2019-06-04 | 2019-11-05 | 浙江省水利河口研究院 | A kind of calculation method of tidal bore flow velocity vertical characteristics |
| CN110555236A (en) * | 2019-07-31 | 2019-12-10 | 浙江省水利河口研究院 | large-size fixed marine platform wave vertical force calculation method |
| CN111553069A (en) * | 2020-04-23 | 2020-08-18 | 河海大学 | Method for measuring balance degree of terrain profile under wave action |
| CN111723461A (en) * | 2020-05-13 | 2020-09-29 | 水利部交通运输部国家能源局南京水利科学研究院 | Wind stress coefficient expression method and system comprehensively considering influences of wind speed, blowing range and water depth |
| CN111783345A (en) * | 2020-07-13 | 2020-10-16 | 中国水利水电科学研究院 | Water ice sand coupling numerical simulation method based on river ice dynamics |
| CN112632792A (en) * | 2020-12-30 | 2021-04-09 | 中交第一航务工程局有限公司 | Shore-approaching wave reanalysis simulation system based on measured data and numerical simulation |
| CN112749520A (en) * | 2021-01-20 | 2021-05-04 | 中国科学院南京地理与湖泊研究所 | Three-dimensional hydrodynamic force numerical model modeling method for large shallow lake |
| CN112784502A (en) * | 2020-12-18 | 2021-05-11 | 清华大学 | Hydrologic-hydraulic dynamic bidirectional coupling flood prediction method |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2008110909A2 (en) * | 2007-03-15 | 2008-09-18 | Universita' Degli Studi Di Palermo | Measurement method for measuring the peak flow discharges and the channel bed roughness in a natural or artificial water course and related system |
| CN105205200A (en) * | 2015-04-27 | 2015-12-30 | 河海大学 | Verification method for silt coast wave-induced current numerical simulation |
| CN107330196A (en) * | 2017-07-03 | 2017-11-07 | 中国地质大学（北京） | It is a kind of to predict the method that the small yardstick coastline configuration of chiltern loke shore is developed |
| CN107895059A (en) * | 2017-09-18 | 2018-04-10 | 水利部交通运输部国家能源局南京水利科学研究院 | A kind of silt coast high concentrtion sea area islands and reefs promote silt engineering simulation method |
| CN108256137A (en) * | 2017-09-18 | 2018-07-06 | 水利部交通运输部国家能源局南京水利科学研究院 | A kind of Macro-tidal estuary gulf man-made island operation area basin sedimentation in channel analogy method |
-
2018
- 2018-08-19 CN CN201810944700.4A patent/CN109582996A/en active Pending
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2008110909A2 (en) * | 2007-03-15 | 2008-09-18 | Universita' Degli Studi Di Palermo | Measurement method for measuring the peak flow discharges and the channel bed roughness in a natural or artificial water course and related system |
| CN105205200A (en) * | 2015-04-27 | 2015-12-30 | 河海大学 | Verification method for silt coast wave-induced current numerical simulation |
| CN107330196A (en) * | 2017-07-03 | 2017-11-07 | 中国地质大学（北京） | It is a kind of to predict the method that the small yardstick coastline configuration of chiltern loke shore is developed |
| CN107895059A (en) * | 2017-09-18 | 2018-04-10 | 水利部交通运输部国家能源局南京水利科学研究院 | A kind of silt coast high concentrtion sea area islands and reefs promote silt engineering simulation method |
| CN108256137A (en) * | 2017-09-18 | 2018-07-06 | 水利部交通运输部国家能源局南京水利科学研究院 | A kind of Macro-tidal estuary gulf man-made island operation area basin sedimentation in channel analogy method |
Non-Patent Citations (5)
| Title |
|---|
| 姬泓宇等: "粤东靖海湾海岸地貌动力演变及其工程影响", 《海洋工程》 * |
| 张海文等: "近岸波、流作用下结构物附近海岸演变的数值模拟", 《海洋学报(中文版)》 * |
| 李国胜等: "大河三角洲河口海岸演化机理模型研究:(I)模式理论与进展", 《地理研究》 * |
| 郭雅琼等: "莺歌海离岸码头接岸方式对岸滩演变影响研究", 《科学技术与工程》 * |
| 陈沈良: "琼州海峡南岸海岸动力地貌研究", 《热带海洋》 * |
Cited By (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110388899B (en) * | 2019-06-04 | 2022-08-09 | 浙江省水利河口研究院 | Special assembly for calculating tidal bore flow velocity vertical distribution |
| CN110412312A (en) * | 2019-06-04 | 2019-11-05 | 浙江省水利河口研究院 | A kind of calculation method of tidal bore flow velocity vertical characteristics |
| CN110412312B (en) * | 2019-06-04 | 2021-05-11 | 浙江省水利河口研究院 | Calculation method for tidal bore flow velocity vertical distribution |
| CN110388899A (en) * | 2019-06-04 | 2019-10-29 | 浙江省水利河口研究院 | The personal module calculated for tidal bore flow velocity vertical characteristics |
| CN110555236A (en) * | 2019-07-31 | 2019-12-10 | 浙江省水利河口研究院 | large-size fixed marine platform wave vertical force calculation method |
| CN111553069A (en) * | 2020-04-23 | 2020-08-18 | 河海大学 | Method for measuring balance degree of terrain profile under wave action |
| CN111553069B (en) * | 2020-04-23 | 2022-11-04 | 河海大学 | Method for measuring balance degree of terrain profile under wave action |
| CN111723461A (en) * | 2020-05-13 | 2020-09-29 | 水利部交通运输部国家能源局南京水利科学研究院 | Wind stress coefficient expression method and system comprehensively considering influences of wind speed, blowing range and water depth |
| CN111723461B (en) * | 2020-05-13 | 2022-08-23 | 水利部交通运输部国家能源局南京水利科学研究院 | Wind stress coefficient expression method and system comprehensively considering influences of wind speed, blowing range and water depth |
| CN111783345A (en) * | 2020-07-13 | 2020-10-16 | 中国水利水电科学研究院 | Water ice sand coupling numerical simulation method based on river ice dynamics |
| CN112784502A (en) * | 2020-12-18 | 2021-05-11 | 清华大学 | Hydrologic-hydraulic dynamic bidirectional coupling flood prediction method |
| CN112784502B (en) * | 2020-12-18 | 2022-09-20 | 清华大学 | Hydrologic-hydraulic dynamic bidirectional coupling flood prediction method |
| CN112632792B (en) * | 2020-12-30 | 2022-03-04 | 中交第一航务工程局有限公司 | Shore-approaching wave reanalysis simulation system based on measured data and numerical simulation |
| CN112632792A (en) * | 2020-12-30 | 2021-04-09 | 中交第一航务工程局有限公司 | Shore-approaching wave reanalysis simulation system based on measured data and numerical simulation |
| CN112749520A (en) * | 2021-01-20 | 2021-05-04 | 中国科学院南京地理与湖泊研究所 | Three-dimensional hydrodynamic force numerical model modeling method for large shallow lake |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109582996A (en) | A kind of coupled simulation method of small scale beach profile and large scale Shoreline changes | |
| CN108256137B (en) | Method for simulating channel desilting of harbor basin in artificial island operation area of strong-tide estuary | |
| Shi et al. | Anthropogenic influences on the tidal prism and water exchanges in Jiaozhou Bay, Qingdao, China | |
| Zhou et al. | Refinements to the EFDC model for predicting the hydro-environmental impacts of a barrage across the Severn Estuary | |
| Hu et al. | Modeling the mass flux budgets of water and suspended sediments for the river network and estuary in the Pearl River Delta, China | |
| Kuang et al. | A two-dimensional morphological model based on next generation circulation solver II: Application to Caofeidian, Bohai Bay, China | |
| CN110135069B (en) | Method and device for acquiring silt characteristics during water delivery of water delivery tunnel and computer equipment | |
| Brown et al. | Methods for medium-term prediction of the net sediment transport by waves and currents in complex coastal regions | |
| CN107256312B (en) | Gulf tide receiving variation calculation method based on tidal current environment | |
| Zhang et al. | Large-scale 3D numerical modelling of flood propagation and sediment transport and operational strategy in the Three Gorges Reservoir, China | |
| Xie et al. | A preliminary study of the turbulence features of the tidal bore in the Qiantang River, China | |
| He et al. | Simulation of sedimentary dynamics in a small-scale estuary: the role of human activities | |
| Zhang et al. | 3D numerical modelling of asynchronous propagation characteristics of flood and sediment peaks in Three Gorges Reservoir | |
| Zhang et al. | Simulation and analysis of Back siltation in a navigation channel using MIKE 21 | |
| Hou et al. | The morphodynamic responses to deposition-promoting projects in island and reef coasts of the Zhoushan Archipelago, China | |
| Zhu et al. | Response of stratification processes to tidal current alteration due to channel narrowing and deepening | |
| Thong et al. | Numerical study of sediment transport in Thu-Bon estuary and coastal areas of Hoi-An City | |
| Lianqing et al. | Numerical simulation and optimal system scheduling on flood diversion and storage in Dongting Basin, China | |
| Li et al. | Numerical simulation of tidal current energy in Yangtze Estuary-Hangzhou Bay, China | |
| Sharbaty | Two Dimensional simulations of seasonal flow patterns in the Gorgan Bay | |
| Ma et al. | Identification and simulation the response of storm-induced coastal erosion in the China Yellow sea | |
| Han et al. | Numerical experiments on stagnation points influenced by the Three Gorges Dam in the Yangtze Estuary | |
| Sharbaty | 3-D simulation flow pattern in the Gorgan Bay in during summer | |
| CN115796077B (en) | Method for configuring irrigation water in tide draining and irrigated area, computer device and storage medium | |
| Zhang et al. | Nested numerical simulation of multi-scale flow characteristics around monopile foundation of offshore wind farm at Weihai coastal water, China |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20190405 |