TW201201035A - Two-dimensional hydraulic sediment transport simulation calculation method and computer program product - Google Patents

Abstract

A two-dimensional hydraulic sediment transport simulation calculation method adopts finite volume method to perform numerical discreteness of a governing equation, which is applied to the unstructured grid to allow hydraulic simulation calculation of irregular river boundary and hydraulic structure in the river. Using finite-volume multi-stage schemes (abbreviated as FMUSTA) for the numerical flux estimation can inhibit numerical oscillation and increase the stability of numerical calculation, and can perform hydraulic simulation of sub/supercritical flow conditions. Using explicit time integration method can perform hydraulic simulation of fixed/variable flow. Using surface gradient method for the treatment of bed slope source term allows dry/wet hydraulic simulation under irregular landform. The FMUSTA calculates inviscid numerical flux to have the relative error reduced, and operation time of the processor is substantially shortened. All in all, it completes simulating two-dimensional hydraulic calculation in the most efficient way, and further simulates the phenomena of sediment transport.

Description

201201035 VI. Description of the Invention: [Technical Field of the Invention] The present invention relates to a numerical simulation calculation method, and in particular to a two-dimensional simulation method for water treatment and sand transport. [Prior Art] The study of the watery phenomena of the shallow water (referred to as the aspect ratio) in the river estuary is of great importance for the protection of land and pollution. The water phenomenon refers to the change of the water flow rate and the water level in the water. It can be simulated in the one-dimensional mode Φ, two-dimensional mode or three-dimensional mode. For alluvial rivers with a small aspect ratio, the depth direction changes much less than the lateral variation, and can be simulated in horizontal two-dimensional mode. In recent years, under the numerical discrete structure of the finite volume method (FVM), using the characteristics of the aerodynamics in the anti-wind conservation calculation, many successful application cases have been used in the open channel hydraulics. Numerical simulation of the flow field. The characteristic upwind conservation algorithm changes the traditional flow rate upwind into appropriate dissipative according to the characteristic upwind, and automatically adjusts the dissipative intensity according to the φ-free oscillation requirement of the solution, so as to prevent the numerical non-physical when calculating the discontinuous solution. Sexual oscillations can increase accuracy. These cases include, for example, "Approximate Riemann solvers in FVM for 2D hydraulic shock wave modeling" published by Zhao et al. in the Journal of Hydraulic Engineering in 1996, using approximate Riemann solvers, including: Osher, Roe, and FVS (flux vector splitting) three kinds of calculations to solve the two-dimensional shallow water wave equation; Brufau and Garcia-Navarro in 2000 at International 201201035

The "2-dimensional dam break flow simulation" published in the journal Journal for Numerical Methods in Fluids, based on the Roe algorithm, uses the slope-limiter approach for the simulation application of the two-dimensional dam test; Macchione and Morelli The ''Practical aspects in comparing shock-capturing schemes for dambreak problems' published in the Journal of Hydraulic Engineering in 2003 compared and explored the difference in precision between many different first-order and second-order algorithms. The results show that for irregular bed slopes And the accuracy of the first-order and second-order calculations is not much different when the shallow-water flow field simulation of the bed resistance is used; Wang et al. used the "2-dimensional free surface flow in branch channels by a finite-volume TVD scheme" in 2003. The flux-limiter approach uses the Roe algorithm to simulate the problem of two-dimensional dam flow passing through the branch channel. The applicant Lai team developed the upwind flux-splitting finite (upstream flux-splitting finite) Volume scheme, referred to as UFF), for two-dimensional shallow water Numerical simulation of the field, and the simulation results are compared and compared with Osher, Roe and HLL calculations. The results show that the numerical precision and numerical performance of the UFF algorithm are slightly better than the other three. The anti-wind conservation algorithm can be successfully applied to solve the shallow water wave equation with the finite volume method, and it has great analytical ability for the shock wave with discontinuous water flow phenomenon, and does not produce non-physical numerical oscillation phenomenon. Most of these calculations are for the simulation of shallow water flow field with regular bed bottom slope. For shallow water flow field model 201201035 with irregular and complex terrain, the source term must be specially treated to avoid flux gradient and The error caused by the non-conservation of the source item. The applicant's Lai team published in December 2005, Taiwan Water Conservancy, Volume 53, No. 4, "The anti-wind flux splitting finite volume algorithm is in the water flow of the Keelung River round mountain section. Simulation"—in the text, based on the Qi Qing calculation, the value of the source term is further extended by the SURface gfadient method (SGM). First, the time-limited volume method is used to control the financial program, and the finite volume equation is obtained. Then, the numerical problem is solved. The I-transformation is a one-dimensional Riemann problem in the vertical direction of the interface between elements and elements. To deal with 'and estimate the numerical flux; supplemented by the use of reading method for the bed source item processing, the water depth data originally used in the control equation & the original water depth data construction, changed to use water level to deal with, successfully improve the upwind flux split type Source processing capability for finite volume calculations. SUMMARY OF THE INVENTION It is therefore an object of the present invention to provide a two-dimensional water transport simulation simulation method for numerical flux estimation using Finite Volume-Volume Multi-Stage Schemes (FMUSTA). The FMUSTA algorithm can reduce the relative error and the calculation efficiency is higher. Another object of the present invention is to provide a computer program product, a system program for storing a two-dimensional water sand transport simulation calculation, and when the computer is loaded into the program and executed, the two-dimensional water transport sand simulation calculation method of the present invention can be completed. . Therefore, the two-dimensional hydraulic sand transport simulation calculation method of the present invention is executed by using a processor of an electronic device to read the code, and comprises the following steps: (a) calculating a first-order convection numerical flux of the adjacent element interface, wherein 201201035 Z FMUSTA algorithm, first converts the local-Villemanian equation used to estimate the non-viscous flux of the vertical boundary of adjacent elements into a local coordinate system, and calculates the two parts: (4) and the correction part. The left element physical vector and the right element physical vector are used as the correction flux by the left element initial vector of the element interface and the right element initial vector as the predictive flux re-matching constant formula to estimate the element interface of the next moment. 'is the numerical fluence of the convection; (b) expands the first-order convection numerical flux of the adjacent element interface to the second-order precision; (c) calculates the diffusion numerical flux of the adjacent element interface; (d) performs the element center (e) Solving the water depth, flow velocity, and sediment concentration of the sediment in the center of the element; (f) Calculating the silt bed load using a riverbed sediment transport formula; (g) Solving the bed Cheng; and (h) outputting the simulation results. Preferably, step (a) is to use a local Lax-Fdeddchs numerical flux function to obtain a predicted flux through one-half of the element interface. Preferably, the step (a) comprises the following substeps: (al) obtaining an initial maximum value of the gravity wave velocity by the water depth and the flow velocity under the local coordinate system; (a2) the physical vector by the local coordinate system And the flux, and the initial maximum value of the gravity wave velocity obtained by the step (al), with a Koran numerical stability 201201035 conditional coefficient, to solve the predicted flux through the element interface; (a3) by the conservation formula and the step (〇The predicted flux is calculated to solve the physical vector at the next moment; U4) using the physical vector at the next moment to obtain the maximum value of the gravity wave velocity; and (a5) is obtained by the step (a3) The next time physical vector and the maximum value of the gravity wave velocity obtained in step U4) are obtained. Preferably, the present invention further includes a pre-processing step φ before step (a), including receiving terrain data, generating a calculation grid within the simulation calculation range, and receiving relevant settings required for the simulation calculation, including at least one of the following: Manning drag coefficient, simulation calculus time or number, numerical simulation time interval, dry and wet processing parameters, sediment-diffusion diffusion related parameters, turbulent motion viscous coefficient, options for four different riverbed sedimentation transport formulas, flow rate, water level And the initial conditions and boundary conditions of the flow rate. Preferably, the relevant settings received by the pre-processing step further comprise receiving inflow and outflow conditions of the water flow and the mud sand, and step (e) also solving the sludge mass concentration of the mud φ sand. Preferably, the present invention further comprises a step (i) performed before or after the step (h), comparing the simulation result with the actual measured data obtained by the additional investigation, and if the coincidence is completed, if not, the Receive the relevant settings required for the simulation calculation. The efficacy of the present invention can be found from the numerical flux function presented above, and the FMUSTA algorithm is more accurate than the Roe algorithm that is often used, and the arithmetic operators, relational operators, operands, and expressions are more accurate. All are less than R〇e 201201035, which saves a lot of calculation time. Therefore, it is quite simple to write a calculation program, so it is worth promoting. The above and other technical contents, features, and advantages of the present invention will be apparent from the following detailed description of the preferred embodiments. The preferred embodiment of the two-dimensional hydraulic sand transport simulation calculation method of the present invention is implemented by a system program using 'C++ syntax, and the development technology used includes Windows Programming and Computer Graphics. And computer Aided Geometry Design, etc., and integrate the hydraulic simulation program as the simulation core. The system program can be pre-recorded in a computer program product. When the processor of an electronic device reads the code and executes it, the two-dimensional water transport sand simulation calculation method of the present invention can be completed. Referring to Fig. 1, the two-dimensional hydraulic sand transport simulation calculation method of the embodiment comprises the following steps: Step S11 - receiving data. The system can import real terrain elevation data such as aerial photographs and other terrain elevation data to make a calculation grid in the subsequent steps, and combine the hydrological sand transport simulation calculation program to complete the numerical simulation work. Step S12—Generate the simulation calculation range The calculation grid inside. The system mainly generates grid according to the imported data, and provides two-dimensional geometric design function, * includes various mode switching, so that the user can add, select, edit attributes and points of points, lines, etc. in the system. Delete and other actions. The system 2010201035 also provides a function for smoothing the structured grid. The system receives the smoothing processing instruction input by the user, receives the iteration number set by the user and smoothing the parameter, and then selects the smoothing parameter. All structured grid blocks regenerate smoothed internal mesh data and smooth it. In addition, the system also provides the function of selecting the peripheral boundary of the grid and the selection of the grid elements, and the re-refinement of the grid. #关网 re-refinement re-creation, for example, to make a more detailed simulation analysis for the scouring condition of the pier Φ. The calculation grid elements in the vicinity of the bridge can be selected to be classified by the connection relationship, and several grid element divisions are divided. Then, each mesh element partition is set in detail, and the setting method is as follows (for each grid element partition): receiving a command of a new control point determined by the user in the divided area, determining a feature The instruction of the line segment (at least 2 control points is required), and then the line segment is uniformly distributed to generate detailed feature points, and the setting of such a feature line segment is completed. • 2. Ask if you want to continue to decide on additional detailed feature lines. If the user chooses to continue, the user will continue to let the user decide the new feature line segment. If the user chooses not to continue, the feature line segment setting of the segmentation area ends. Continue to make detailed settings for the next segmentation area. When the detailed settings of all the grid element partitions are completed, the system starts the internal mesh re-creation operation for each partition, and displays the re-made results as a preview to the user as shown in FIG. Ask the user whether to accept the result of this grid re-creation. If you agree, replace the calculated grid space with the new reworked grid result; 201201035 If you do not agree, restore the original calculated grid space data. Step S13—Receiving relevant settings required for simulation calculation, including Manning resistance coefficient, simulation calculation time or number of times, numerical simulation time interval, dry and wet processing parameters, sediment-diffusion diffusion related parameters, turbulent motion viscosity coefficient, and four types The choice of different riverbed sediment transport formulas, in the case of hydraulic simulation, this step should receive the initial conditions and boundary conditions of flow, water level and flow rate. • 'In the case of sand transport simulation', this step should receive the inflow of water and mud sand. With outflow conditions. In other words, if the water conditions and the inflow/outflow conditions of the silt are input in this step, the water and sand transport simulation results can be obtained at the end; if only the water conditions are input, only the hydrological simulation results are obtained; Conditions and mud sand inflow/outflow conditions, and finally the results of water and sand transport simulation. When the grid is completed and the conditions are set, the simulation flow shown in Figure 3 is executed. Before explaining the actual calculation flow, the following is a brief description of the theoretical derivation part of the present invention. The present invention uses the finite volume method to perform numerical dispersion of the control equation. First, in order to solve the hydrological phenomenon of sediment transport in the shallow water of the river, combined with two-dimensional The shallow water wave equation and the silt suspension mass transfer equation are written as the following vector conservation type: (Formula 1) + j.9 c〇nv - ^Fdiff BG diff Seven dx ~dy~ where, 10 201201035

•/Γ hu " Q- hu ;F = hu2+〇.5gh2 hv Iconv huv hC _ huC G conv hv Ί huv as 2+0.5 zhen

hvC ' 0 ' ' 0 ' P Ά ' G<Mff = P hT. P • 1 υχ^ψ j^d(hQ ;s= 8h{S(>y~sfy) -%~qsd 0

, Q is the prune physics vector; f, g 八 丨 convection (_ sex) flux vector; ', ^ respectively ^ is the ^ direction (viscous) flux vector; s is the source term, 1 contains the diffusion : Suspension - load source ... water depth; ... respectively for the direction of the water depth average '• speed for gravity acceleration; ~=_</such as, \=_3 ~ / do and L, ~ respectively named X' 3; the bottom of the direction Bed slope and friction gradient; 4 is bed height; I, && is horizontal turbulent stress; p is fluid density; 匚 is water depth average silt suspension concentration; A and A are diffusion coefficients; The turbulent shear force produces the flux of the riverbed upwards; "the flux that sinks down to the bed surface for the suspended matter. The friction gradient can be expressed as follows: (Formula 2) S{x=^^u2+f_ . ^ = v In the "annV" 2 + v2 formula, it is the Manning coarse wheel coefficient. In addition, according to Boussinesq's 11 201201035 thirst viscous concept, the flow stress can be written in the form of a component viscous term, then Γ , ^ and L are expressed as follows: λ

B(hu) η/. λ/.-,Τ-ν: ?-~ί~ ^(hv). ~βΓ]'

(Formula 3) where ' is (4) the viscous coefficient of motion, which can be estimated by the empirical formula ^"(4); where '-_ coefficient, using 〇.4... is the shear velocity; Tbx=PghSfiM τ is χ" Bottom bed shear stress. Diffusion coefficient, can be public 戎 n _ n „ . ^ A type of estimation, where σ is turbulent

The Prandtl-Schmidt coefficient needs to be controlled between 〇 5 and j 〇. Next, the numerical dispersion of the governing equation is performed using the finite volume method of the element center type. It is assumed that each variable of the physical vector Q is a physical variable that controls the geometric center point (centroid) of the volume area. For any control body Ω in the simulation calculus domain, the control equation is integrated, and the divergence theorem is used to differentiate the area of the physical flux into a line integral along its perimeter boundary:

3Q dt dW +

ilSdw (style 4)

Where dW is the microvariable of the control body area · E - (F, G) - (Fe(jnv - , GeDnv - Gjjjf ) is the flux vector of the text, ^>plane; the buckle is around the body Ω π is a unit vector that is positive to the outside of the boundary 3Q. dZ is the microvariable of the length of the control body. According to the assumption of the variable of the element center, 12 201201035 and the rotation invariance of the governing equation, the equation 4玎 is further discretized into a finite volume. The basic algebraic equation of operation: ^ + TEWF(Q)Lffl= S (Formula 5) where 'A is the area of the control body; Μ is the total number of perimeter boundaries of the control body; L is the length of the wth side of the control body; The angle between the vector η and the X axis (calculated from the X axis in the counterclockwise direction); the physical vector Q in the X '^^^ system is transformed by the coordinate rotation to obtain the conversion physics vector under the 7K system in the new block of r, y Q, where the new coordinate system is the vertical boundary ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The new seat is taken, the dead squat, the flow velocity in the 歹 direction; F along) = Fe()nv(7-F For the vertical boundary flux after transposition; Τ(6) and τπγ are the transposed matrix and the inverse transposed matrix respectively 疋 as follows: '1 0 0 0' 'l 0 0 0 — Τ(6>) = 0 COS0 sin0 0 0 COS0 -sin0 0 0 ~sin0 cos^ 0 ,i (β)- 0 sin0 COS0 0 _0 0 0 1 0 0 0 1 (Formula 6) On the left side of Equation 5, the first term on the left side of the equal sign represents the time that has not been separated. Item, etc. The two terms represent the spatial discrete term 'the right side of the equal sign represents the non-homogeneous source term 13 201201035. This embodiment uses the method of predicting and correcting the Muscx (monotonic upstream schemes for conservation laws) in the subsequent step S6, and obtains the long solution The space and time are the second-order precision calculations: , wr drink - is the name of the heart (_ + side (style 7a) A ·,; Σ τ (0) " ^ ν (〇) ^ + Σ WFdiff (Q) r +s _m=lm=t n+1/2 (Equation 7b) where y· represents the discretized spatial index of the λ, y direction, respectively; a, is the area of the element α)); ql represents the moment The element 0·, is the physical vector. Equation 7: First, the physical vector q under the ;c,y coordinate system is transformed into the physical vector $ under the new coordinate system via the transposed matrix, and F(G) is calculated by ^; then the element length and reversal are multiplied. The matrix is changed to the vertical flux under the original coordinate system, and the sum of the vertical fluxes of the sides of the element, plus the source term processing, gives the final numerical solution. The spatial discrete term in Equation 7 includes the convective numerical flux (strand) and the diffusion numerical flux Fditf (5). In step S2, the estimation of the order convection numerical flux of the adjacent element interface is included. Sub-step S21 § 25 is shown as for the diffusion numerical flux aspect 'using the commonly used central difference solution, which will be performed in step S4. The Fdiff (G) through the element interface is: 14 201201035

Fdi«(Q) = 2/4 0 3(/zWj)

9(3⁄4 & dp M 〇RD ^(hQ • X~W (Equation 8) is where the gradient derivative term of the velocity and sediment concentration can be expressed (Equation 9) (Formula 10) jOO _[(^), >ιι7· ~ ](y,.,.|/2,j+l ~ ^,+1/2.7-1) ~[(·),-.-1/2,j+i ~ ~p, ί ( ^+1./ ~ (*)|,;](-^|+1/2./+1 -文丨+丨...-丨)-[(·),+丨/2,j+l - (·) ι·+ι/2,)-丨] (also Hl.j -',) (^•+1,;· - yij)(xMi2,j+\~^i+mj-i)-(yMi2j +i - yi+in,ji)^Mj~^ij) where the physical quantity in the brackets can be; 1Μί, /jVy or . Back to step S2, in the estimation of the convective numerical flux uu), The physical vector values of each element in the computational domain may have different values, and thus may be discontinuous when the interface between elements and elements. Because of this discontinuity, when estimating the non-viscous flux at the boundary of each element, this problem can be transformed into a local one-dimensional Riemann problem in the vertical direction of the element-to-element interface: M" column U] 93c < ο > ο (Expression 11) where χ is the origin of the element in a certain boundary - and the vertical outward is a positive coordinate axis (as shown in Figure 5); F_ (received) is a vertical from the i-axis origin The outward flux is called the = non-viscous flux vector. The initial condition of Equation 11 is defined as: In the absence of <〇, the temple is 仏, where 屯 and 仏 are the physical vectors to the left and right of the element 兀 interface LR, the left element is the element being calculated, the right Elements are adjacent elements. A variable (pure 15 201201035 quantity or vector) whose subscript L, and can represent the left element L and the right element and the variable # value. Suppose $ is in? = G is a known value, so the Riemann problem of the initial value can be solved to find the non-viscous flux of the vertical boundary at r=0+. The present invention is characterized in that a Finite volume multi-stage scheme is applied to (Formula 11) and converted into the following partial coordinate system'-- and step S21 of the present embodiment calculates the following (Formula 12 S+) &wS)i dt dd :〇, Q(rf,〇) = 〇T = QL d<0 [q:0) = qr J>0 (Expression 12) where 7 is a spatial variable relative to X; Representing the time variable relative to £, the subscripts and i of the variables represent the initial vector of the left element 〇 and the right element i of the element interface 1/2, respectively, relative to the elements L and R: (f) and Q, respectively. Shown. The calculation is divided into two parts: prediction and correction. Step S21 and SU of this embodiment perform the prediction part, which uses the lGeal Lax Fdeddehs numerical flux function formula to obtain the predicted flux through the element interface 1/2: (Expression 13) where CFL is a kolan numerical stability condition Coefficient (c〇urant_Friedrichs-Lewy stability coefficient); ^ t Λ ^ ^ it ^ ^ The initial maximum value can be expressed as follows in step S21, and then in step S22, the pass element 13 is taken back to find the pass element The predicted flux of the interface. Max (〇) 〇 (Expression 14) Next, in step S23, the physical vector of the next-time is obtained by the conservation formula and the initial prediction flux 16 201201035: (Equation 15a) (Equation 15b)

Wherein 'private' and ^0 respectively represent the left element 〇 of the element interface 1/2 and the (1)th time physical vector of the right element 1. Finally, the flux estimation of the correction part is performed in steps S24 and S25: In the formula, the maximum value of the (1) time of the gravity wave velocity can be displayed as

(Equation 17) The corrected step flux in Lu 16 is the convective numerical flux.

LuSWr). Next, the first-order numerical flux is expanded to the second-order precision in step S3. The calculation of the second-order numerical flux is based on the MUSCL method. It is assumed that the physical variables are linearly distributed within the control body element, and the weighted linear interpolation is performed according to the slope of the variable (variable difference) to obtain the physical quantity adjacent to both sides of the interface of the element, which is called the variable. Construction. In order to be consistent with the finite volume discretization, the conservation variables can be used to construct the left and right physical quantities of the element boundary surface 0 + 1/2,)) as follows: 17 201201035 1/2,; (Equation 18a) (Equation 18b) See (8)= W) (Q, _+u where 觅丨/2,; and 丨/2,; represent the left and right physical quantities of the element boundary surface (ι·+1/2, y·); Δυ =4; (Δ,.+1/2;,Δ,·_Ι/2;) is an element (/, is a two-parameter slope limit subfunction. In order to avoid excessive slope during interpolation, the numerical solution may oscillate, so it needs to be Limit slope: ij [sgn(A-+l^2 )+sgn(A(i_^2;· )]τ |Δ<'+1/2.;ΊΔ«'-1/2,;| |^i +l/2,)| + Ai-l/2.J + f (Equation 19) where Δ.+1~=Qi+1J -CM port ΔΜ~_ =Q,. Wide QMJ is the variable difference (slope) Sgn represents the sign function, f is a small positive number. Then, the left and right physical f 〇u of the element side, 丨 (ί+ι/2,) can be entered into the FMUSTA numerical flux calculation formula. The numerical flux has a second-order accuracy. The convective numerical flux is calculated in steps S2 and S3, the diffusion value is calculated in step S4, and After the recursive calculation for the adjacent element surface treatment, the process proceeds to step S5 to perform the processing of the source item of the bottom slope. 〇 π This embodiment uses the water level minus the bottom elevation at the element boundary surface to wait for the interface left and right. Water depth on both sides: 18 201201035 hLn.j

VI ϊ+1/2,7 — ^fei+l/2,> J ^1+1/2,> = 7i+l/2,; ^ Zbi (Formula 20) where 仏~ and heart are respectively The boundary surface of the representative element "+1/2, the water level on the left and right sides of the knife" can be obtained by the construction of (10); the boundary surface of the element is G+1/2, the elevation of the bottom of the force; 吣~ and 仏... represent the boundary surface of the element The depth of the left and right sides of the water ^ final 'can directly use the divergence theorem to perform the numerical dispersion of the two terms: ' (Formula 2la) 'J h-w2m; c4)(yi-y3) + —~ Ί. ii, j (Eq. 21c) Four nodes (in the inverse time formula, 'subscripts 1 to 4 respectively represent the element α)) clock direction number). After the calculation is completed, the element can be solved in step S6. In the middle, the flow rate tt, v and the sediment concentration of the mud sand. T love in the step... will carry out the bed bed _ calculation. - This is the solution caused by the continuous equation and the suspended carrier transport, the equation can be It is expressed as follows: /7 Bed deformation around the structure in Sichuan (Equation 22) 19 201201035 Where, the meaning is the porosity of the sedimentation; ^ The density of the sedimentation rate of the riverbed in the direction of the sedimentation rate of the riverbed sedimentation rate is as follows (1) Shamov Formula: In step S7, the water depth and flow obtained by (4) s6 are substituted into the following four types of riverbed sediment transport formulas, and the sedimentation rate is calculated by the calculation of the sedimentation rate and the sedimentation rate. The calculation formula is respectively converged (Equation 23), where % is the transport rate of riverbed sedimentation; & is the empirical coefficient; it is the median granule of sedimentation; 7 is the average flow velocity; R is the starting velocity of muddy sand, which can be expressed as:

Ur = 1.34 rs-r 'γ noisy φ. 143 (Expression 24) where h and y are the specific gravity of mud and water, respectively. (2) Meyer-Peter & Muller (MPM) formula qb = 8 (-) 05 [j / w / -0.047 (^ - /) D] (21zZ)

Y (Equation 25) (3) van Rijn Formula: qb=〇.〇53^s/Y-\)gDl5^t where ϋ is the motion viscous coefficient; £>. and Γ* are defined as follows (Equation 26) υ pghSf ~0.05(ys-y)p 0.05(~y)D ] (Equation 27) 20 201201035 (4) Young's formula:

wD U log ¢6 = 5.435 -0.286 log(-)0.457 log(-) vw + [1.799-0.409log(-)-0.3141og(^)]log[^--r VWW vv JV ^o) where w is Settlement velocity; c/· is the friction speed; 匕 is the critical starting speed.

By defining Re* = ί/*Ζ)/υ as the Reynolds number, the critical starting speed can be expressed as follows:

Uc 2.05 (70 幺 Re.) (Equation 29a)

U

2.5 log(Ret )-0.06 + 0.66 (1.2 < Re, < 70) (Equation 29b) In addition, the sedimentation velocity can be expressed as: VV = 2i6 (Z7Z) gV (D < 01mm) (Equation 30a) 1 (13.95 I) 2+1·〇9〇^Ι 7

) sand-13.95 upper D (form 30b) 21 201201035 (formula 30c) The numerical dispersion is also based on the finite volume method: (Equation 31) Next, the divergence theorem is used:

From & : 1 A M (Formula 32), “The transport rate of the riverbed through the (10) boundary element can be selected from different sand transport formulas. After calculating the sediment carrying rate of the riverbed in step S7, the step of carrying the sediment carrying rate of the riverbed and the silt suspension loading (that is, the mud calculated in step %: the concentration of suspended matter C) is brought back to the bottom of the solution. Bed elevation changes. Finally, 'Yes 6 is not completed, if it is not completed, then step S2~S8 h 1 will be completed with 70% of the elements, and then confirm whether the number of executions reaches the preset has not reached the preset reading. The element is executed repeatedly... /8'm is counted up to the preset number of times, and the simulation result is obtained. Refer to circle 1. When the above steps are performed to complete the simulation, the output is output for the simulation results. In this embodiment, the simulation result is mainly presented in a graphical manner by snoring the drawing software, and numerical solutions such as water level, water depth, current speed, mud 9 suspension concentration and bed shape change are followed, and step S15 is less. The simulation results are compared with the measured data. 22 0 201201035 Verify the accuracy and rationality of the system; the measured data refers to additional on-site measurement or collection of existing hydraulic and silt data from previous literature. If the simulation result matches the measured data, the process ends; if it does not match, return to step S12, regenerate the grid, reset the data, or simply reset the data and perform the simulation calculation. The foregoing step S15 is not essential, and the flow may be ended after the completion of step S14. As described above, the present invention uses the finite volume method to perform numerical dispersion of the governing equations, and is applicable to non-structural grids, which facilitates irregular river boundaries and contains hydraulic structures in rivers, such as sand control dams (Fig. 7), Dam and river transport (10) 8), fish pass (Fig. 7 and Fig. 8), pier (Fig. 9), spur (Fig. 1), etc., and the simulation of water and sediment transport; the total variable reduction method can suppress numerical oscillation 'And increase the stability of the numerical calculation, can be applied to the mixed flow conditions of the river flow condition super-subcritical flow, and the discontinuous hydrological phenomenon of the transcritical flow (such as the oblique shock wave, the dam break wave); Time integration method, so the hydraulic simulation of constant/variable flow can be carried out; the numerical flux estimation is carried out by using fmusta algorithm, which is beneficial to dry/wet water simulation under irregular terrain; most importantly, the invention adopts FMUSTA calculates the non-viscous numerical flux. Compared with Roe and other methods, the relative error is reduced, that is, the precision is increased, and the arithmetic operator, relational operator, operand and arithmetic expression of FMUSTA are less than Roe. , Time processor to perform operations also is shortened, the whole 'may use the most effective manner a two-dimensional molded helmet water treatment and sand transport calculations, it can really achieve the object of the present invention. The above is only a preferred embodiment of the present invention, and is not intended to limit the scope of the present invention, that is, the simple equivalent change of the scope of the patent 23 201201035 and the description of the invention. And Senna, are still within the scope of the invention patent. BRIEF DESCRIPTION OF THE DRAWINGS FIG. 1 is a flow chart of a preferred embodiment of a two-dimensional hydraulic sand transport simulation calculation method of the present invention; FIG. 2 shows an internal grid weight for each partition. FIG. 3 is a flow chart of the main simulation calculation program in FIG. 1; FIG. 4 is a detailed sub-step flow chart of step S2 in FIG. 3; FIG. 5 is a schematic diagram of a numerical discrete architecture of the finite volume method; Figure 6 is a schematic diagram of the basic structure of the FMUSTA algorithm; Figure 7 is a schematic diagram of the simulation results of the flow field of the Beishuixi sand control dam containing the fish pass. Figure 8 is a schematic diagram showing the simulation results of the flow rate equivalent distribution of the fish channel system. Schematic diagram of the simulation results of the flow field near the bridge pier in the round mountain section of the Keelung River; and Fig. 1 is a schematic diagram showing the simulation results of the bed bed scouring and silting evolution of the Dingba River section in the Shuiweiwan curve of the Keelung River. 201201035 [Explanation of main component symbols] S11-S15. Procedure S21-S25·Step S 2 - S 8 · · · · Step

25

Claims (1)

201201035 VII. Scope of Application for Patent·· 1. A two-dimensional simulation method for water transport sand transport, which is executed by a processor reading an electronic device, including the following. (a) Calculating the eccentricity of the adjacent element interface Numerical flux, in which the local one-dimensional Riemannian transformation used to estimate the non-viscous flux of the vertical boundary of adjacent elements is first converted into a local coordinate system, and the calculation is divided into two parts: prediction and correction. The prediction part is obtained first. The left-prime initial vector and the right-element initial vector of the element interface are used as the predictive flux, and then the ==: the left element of the element interface of the time-time is calculated: the number of streams: flux physics The vector is used as the correction flux, that is, the -order pair precision (1) expands the first-order convection numerical flux of the adjacent element interface to the second order (1) calculates the diffusion numerical flux of the adjacent element interface; (d) performs the source term calculation of the element center Moisture (1) Cut out the elements (4) Water depth, flow rate and silt suspended matter = Shen: The transport formula calculates the silt bed load; (h) The simulation results are taken out. I π: 利:: A sand transport simulation calculation as described in item 1. The Her-value flux is sufficient to use the bureaux-Freit column flux. (4) to obtain a prediction through one-half of the element interface. 201201035 3. According to the method of claim 2, the two-dimensional hydraulic sand transport simulation calculation method, wherein the step (a) comprises the following sub-steps: a 1) obtaining the initial maximum value of the gravity wave velocity by the water depth and the flow velocity under the local coordinate system; (a2) the gravity wave obtained by the physical vector and flux under the local coordinate system and the step (al) The initial maximum value of the wave velocity is matched with a Coulomb numerical stability condition coefficient to solve the predicted flux through the element interface; (a3) Solving the physics at the next moment by the conservation formula and the predicted φ flux obtained in the step (a2) a vector; (a4) using the physical vector of the next moment to obtain a maximum value of the gravity wave velocity; and (a5) obtaining the next-time physical vector obtained by the step (a3) and the step (a4) The maximum value of the gravity wave velocity is obtained by obtaining the corrected flux. The two-dimensional hydraulic sand transport simulation calculation method according to any one of the claims of the invention, further comprising the pre-processing step before the step (a), comprising receiving the topographic data and generating the simulation calculation range Calculate the grid and the relevant settings required to receive the simulation calculus, including at least one of the following: Manning drag coefficient, simulation calculus time or number, numerical simulation time interval 'dry fish processing parameters, mud sand suspended matter diffusion related parameters, turbulent motion Viscosity coefficient, options for four different riverbed sediment transport formulas, initial conditions and boundary conditions for flow, water level and flow rate. The two-dimensional hydraulic sand transport simulation calculation method described in item 4 of the scope of the patent application 'in which the relevant settings received by the pre-processing step further includes 27 201201035 receiving water flow and mud sand inflow and outflow conditions, and the steps ( 〇 also solves the sediment concentration of the muddy sand and the steps (1) and steps (the elevation of the river bed caused by the muddy sand transfer. Transfer 6. The two-dimensional hydraulic sand transport simulation calculation according to the fourth paragraph of the patent scope The method 'also includes a step performed before or after step (h) 'Comparing the simulation results with the data obtained by the additional investigations. 'If the agreement is the same, the end is 'if it does not match', then the relevant requirements for the simulation calculation are re-received. 7. 7. Computer program product, system program for storing 2D water sand transport simulation calculus * After loading the program and executing it, you can complete the method described in item 1 to 6 of the patent application scope. . 28