CN101956382A - Method for simulating flood routing by solving Saint-Venant equations through channel storage increment relationship - Google Patents

Abstract

The invention relates to a method for simulating flood routing by solving Saint-Venant equations through a channel storage increment relationship, which comprises the following steps of: simplifying Saint-Venant basic simultaneous differential equations into a water balance equation and dynamic equations in delta t time interval, and respectively calculating the channel storage increment of a reach according to the water balance equation and the dynamic equations, namely calculating the channel storage increment according to the water balance equation and calculating the channel storage increment according to two dynamic equations at the beginning of and at the end of the time interval. The channel storage increments of the reach calculated by the two methods are equal, so a new and reliable method is provided for simulating the flood routing by solving the Saint-Venant equations by a hydraulic method and plays an important role in the aspects of reasonable determination of flood control project scale, flood forecasting, scientific flood dispatching and the like in river channels, channels, flood diversion and storage areas and reservoirs.

Description

Hold the method that increment relation is found the solution St.Venant equation group simulation advance of freshet with groove

Technical field

The present invention relates to the flood control and disaster reduction science and technology field, particularly the groove according to water balance equation and kinetic equation formula holds the method that increment relation is found the solution St.Venant equation group simulation advance of freshet.

Background technology

China is the frequent country that takes place of a flood, the population in the whole nation 1/2, and 1/3 arable land and main metropolis all are in below the river flood position, and the regional industrial or agricultural output value that threatened by flood accounts for the whole nation 2/3.To alleviate flood damage and utilize the purpose of water resource in order to reach for national economy and social development service, set up the Mathematical Modeling of simulation advance of freshet, in the river course, channel, branch flood storage areas and reservoir determine that rationally aspects such as flood protection works scale, flood forecasting, science scheduling flood are bringing into play important effect.

Describe the single order Quasilinear Hyperbolic type partial differential equation of the rivers and canals unsteady flow characteristics of motion, at first set up by French scholar's St.Venant in 1871.In the ordinary course of things, can't obtain its general analytic solutions.At present, can only use approximate computational methods to find the solution.Approximate computational methods are broadly divided into two classes:

1, hydraulic method:

The direct differential method is exactly a kind of method wherein.The direct differential method is to replace partial derivative with partial difference quotient, and fundamental differential is dispersed turns to difference equation, asks the method that each node numerical approximation is separated on independent variable territory x～t plane difference gridding.Partial derivative in the fundamental differential can be that difference scheme replaces with multi-form partial difference quotient, thereby obtains different difference equations.But from present technological means, achievement in research generally is used for the planning and designing stage, does not adopt as yet in the River Basin Hydrology forecasting model.

2, hydrologic method:

Hydrologic method such as famous " Muskingun method " are the fundamental differential group to be reduced to the water balance equation respectively and to be approximately to hold let out equation, close to separate the water balance equation and hold with difference form and let out equation.Shortcoming is: 1. generalization is thicker relatively, is difficult to complicated current situation and principal character in the labor river course; 2. be difficult to provide the spatial and temporal distributions and the evolution process of dynamic analog flood in the section of water level in the section.

Summary of the invention

Of the present inventionly hold the method that increment relation is found the solution St.Venant equation group simulation advance of freshet with groove, be with St.Venant fundamental differential group after the Δ t period is reduced to water balance equation and kinetic equation formula, can calculate the section groove respectively according to water balance equation and kinetic equation formula and hold increment: the one, can calculate groove according to the water balance equation and hold increment, the 2nd, according at the beginning of the period, two kinetic equation formulas at end also can calculate groove and hold increment.Because the section groove of two kinds of method calculating holds increment and should equate, therefore, this invention is just found the solution St.Venant equation group simulation advance of freshet for hydraulic method a kind of new, reliable method is provided.

Of the present inventionly hold the method that increment relation is found the solution St.Venant equation group simulation advance of freshet with groove, platform branch flood storage areas of the Du family feasibility study report of having carried out in Hubei Province's water conservancy and hydropower surveying and designing institute in late November, 2009 is repaiied in the volume work and has been obtained Preliminary Applications, calculate by the simulation advance of freshet, found that the Zhi Hongdao project (containing the Jiangkou control sluice) that there is no need to build former report, if implement, can reduce investment outlay surplus in the of 200,000,000 yuan.

Of the present inventionly hold the method that increment relation is found the solution St.Venant equation group simulation advance of freshet, belong to the hydraulic method category with groove.This method both had similar place to hydrologic method with existing hydraulic method, and different places is also arranged.

The place similar to hydrologic method is all to have related to the section groove to hold content; Hydrologic method that different is as " Muskingun method " be with the kinetic equation formula approximate become to hold separate with water balance equation connection after letting out equation, and this method is directly the water balance equation to be joined with the kinetic equation formula to separate.

Compare with existing hydraulic method, something in common is all to have related to need the river topography data and according to hydrological data calibration section roughness; Different is, this method with fundamental differential the Δ t period be reduced to the water balance equation and at the beginning of the period, join after two kinetic equation formulas at end and separate.

Concrete grammar is as follows:

Now St.Venant fundamental differential group is reduced to water balance equation (1), kinetic equation formula (2) and kinetic equation formula (3) in the Δ t period:

$\mathrm{\Δ}{Z}_1=\frac{{\stackrel{\‾}{Q_1}}^2}{{\stackrel{\‾}{K_1}}^2}L---\left(2\right)$

$\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="Z"\; filenames="179555DEST\_PATH\_HSB00000208631500031.png"\; state="\{\{state\}\}">Z</figure-callout>}}_2=\frac{{\stackrel{\&OverBar;}{{\mathrm{<figure-callout\; id="2"\; label="Q"\; filenames="179555DEST\_PATH\_HSB00000208631500031.png"\; state="\{\{state\}\}">Q</figure-callout>}}_2}}^2}{{\stackrel{\&OverBar;}{{\mathrm{<figure-callout\; id="2"\; label="K"\; filenames="179555DEST\_PATH\_HSB00000208631500031.png"\; state="\{\{state\}\}">K</figure-callout>}}_2}}^2}L---\left(3\right)$

More than various in: L is section degree long (m); Δ t is period (s); Z is water level (m); Q is flow (m ³/ s); Q is interval flow (m ³/ s); W is section groove water retention capacity (m ³); Δ Z calculates the water-head (m) of section for the upstream and downstream, section;

Be the section average discharge, Δ W=(W2-W1);

Be section average discharge modulus; Each foot-measuring annotate 1,2 be respectively the period at the beginning of, moment at end.

To join the method for separating formula (1), (2), (3) in order illustrating, to be depicted as the figure that the section groove holds increment Delta W (water balance) and Δ W (power) according to above-mentioned equation group.

As shown in Figure 2, according to the figure that water balance equation (1) is drawn, groove holds the difference that increment Delta W (water balance) is the total inbound traffics in section and total outflow in period Δ t, promptly by the process Q that becomes a mandarin _{Last 1}～Q _{Last 2}With go out to flow through journey Q _{Following 1}～Q _{Following 2}The water body that is surrounded.

Draw according to kinetic equation formula (2), (3) as shown in Figure 3, it is at moment t that groove holds increment Delta W (power) ₁Water surface curve Z _{Last 1}～Z _{Following 1}With moment t ₂Water surface curve Z _{Last 2}～Z _{Following 2}The water body that is surrounded promptly calculates section respectively at moment t by upstream and downstream ₁And t ₂Water level convert the difference of flow area to and multiply by the volume that section length draws after average.

When known conditions is Q _{Last 1}, Z _{Last 1}, Q _{Following 1}, Z _{Following 1}, Q _{Last 2}And Q _{Following 2}(or Z _{Following 2}) time, can find the solution Z according to the relation that groove holds increment Delta W _{Last 2}And Q _{Following 2}(or Z _{Following 2}):

According to equation (1), as supposition Q _{Following 2}(or Z _{Following 2}) value, just can calculate corresponding section groove and hold increment Delta W (water balance).

According to equation (2), (3), continue to utilize the Q of above-mentioned supposition _{Following 2}(or Z _{Following 2}) value, just can calculate the section groove and hold increment Delta W (power).

Press above-mentioned relation by computer program and regulate Q _{Following 2}(or Z _{Following 2}) value, it is equal substantially that the section groove that dual mode is calculated holds increment Delta W, the Z of this moment _{Last 2}And Q _{Following 2}(or Z _{Following 2}) be solution of equations.

When calculating section number is n, then need to increase n-2 section flow rate calculation relational expression solving equation group.

As the i=1 that is numbered with the upstream initial section, the section numbering toward the downstream then is i=2, i=3 successively ... i=n-2, i=n-1, i=n.

When section numbering 1＜i＜n, according to water balance equation (1), hold increment as the flow of section and section groove and represent with two-dimensional array, then between section i-1 and the section i, the water balance relational expression between section i and the section i+1 is respectively:

Q(i-1，1)+Q(i-1，2)-Q(i，1)-Q(i，2)＝2ΔW(i-1，i)/Δt---------------(4)

Q(i，1)+Q(i，2)-Q(i+1，1)-Q(i+1，2)＝2ΔW(i，i+1)/Δt---------------(5)

Join solve an equation (4), (5):

Q(i，2)＝(ΔW(i，i+1)-ΔW(i-1，i))/Δt+(Q(i-1，2)+Q(i+1，2))/2+(Q(i-1，1)+Q(i+1，1)-2Q(i，1))/2------------------------------(6)

Also can be similar to the accepted way of doing sth forms such as (7) of writing:

Q(i，2)＝(ΔW(i，i+1)-ΔW(i-1，i))/2/Δt+(Q(i-1，2)+Q(i+1，2))/2------(7)

In the formula: first in the two-dimensional array flow Q bracket is the section numbering, second at the beginning of the period Δ t, last code name (just be 1, end be 2); First of holding in the increment Delta W bracket of two-dimensional array groove is upstream, section section numbering, and second be that downstream, section section is numbered.

Therefore, the full section accumulated slot that calculates by the water balance equation holds increment ∑ Δ W (water balance) and is:

∑ Δ W (water balance)=Δ t (Q (1,1)+Q (1,2)-Q (n, 1)-Q (n, 2))/2--------------(8)

The full section accumulated slot that calculates by the power equation holds increment ∑ Δ W (power) and is:

Description of drawings

Fig. 1 for of the present invention with groove hold increment relation find the solution the method supposition roughness coefficien of St.Venant equation group simulation advance of freshet be 0.031 and 0.035 o'clock the nineteen eighty-three flood tame platform lock of shutting out calculate and measured water level graph comparison diagram;

Fig. 2 holds increment relation with groove and finds the solution the figure that the method for St.Venant equation group simulation advance of freshet is drawn according to the water balance equation for of the present invention;

Fig. 3 holds increment relation with groove and finds the solution the figure that the method for St.Venant equation group simulation advance of freshet is drawn according to kinetic equation formula (2), (3) for of the present invention.

The specific embodiment

Below hold increment relation with groove and find the solution the specific embodiment of the method for St.Venant equation group simulation advance of freshet and be described in further detail of the present invention in conjunction with Fig. 1, but the present invention is not limited to following examples.

Embodiment one:

Hold increment relation with groove and find the solution nineteen eighty-three flood that the method for St.Venant equation group simulation advance of freshet the draws tame platform lock of shutting out and calculate and measured water level graph comparison diagram according to Fig. 1 is described, analyze the tame platform lock of Du to the roughness coefficien of week during group's section nineteen eighty-three flood

The tame platform lock of shutting out is about 21km to week group section (being big vast road), and the start-stop place is respectively tame Tai Zha of Du and week group, and big vast road is the passage that guiding Han River flood enters branch flood storage areas.

Adopt and surveyed 28 in Hong Dao cross section in 1998, measure section and be numbered 1～28, certainty of measurement is levels such as IV, and the section elevation is a Huanghai Sea basal plane.

The observed stage of shutting out under the tame platform lock and helping in week is for freezing Wusong basal plane, and conversion relation is: freeze Wusong elevation=Huanghai Sea elevation+1.85m.The flood data time be 8 days 7 of October nineteen eighty-three up to 9 o'clock on the 15th, the water-level observation period is per hour 1 time, amounts to 171 groups of data.The discharge process of section under the tame platform lock of shutting out calculates by flow formula according to the lock upstream and downstream water level of observation and gets.

The shut out technology path of tame platform branch flood storage areas flood road roughness coefficien of analytical calculation is:

At first should suppose the roughness coefficien of section, be 0.035 as the roughness coefficien of supposing in Table 1;

Hold the flow that principle that increment ∑ Δ W (water balance) and ∑ Δ W (power) should equate is substantially supposed big vast road downstream end Zhou Bangzhan according to the section accumulated slot, the flow that is numbered 28 (week groups) as supposition in 11 o'clock on the 13rd October of nineteen eighty-three section in Table 1 is 2238.14 (m ³/ s) time, iterative computation accumulative total section groove holds increment ∑ Δ W (water balance)/Δ t and ∑ Δ W (power)/Δ t is-351.39 (m ³/ s), it is-0.49 (m that the groove of section 27 to 28 holds increment Delta W (power)/Δ t ³/ s).

Section numbering less than 28 (week groups) greater than the flow of the calculating section of 1 (under the tame platform lock of shutting out) according to formula (7) iterative computation, to number 27 flow be 2238.32 (m as 11 o'clock on the 13rd October of nineteen eighty-three in Table 1 calculating section by formula (7) ³/ s);

Section numbering is finished according to formula (3) iterative computation less than the water level of 28 calculating section, and calculating section by formula (3) in 11 o'clock on the 13rd October of nineteen eighty-three in Table 1, to number 27 water level be 28.12m;

After finishing aforementioned calculation, if under the situation of supposition section roughness coefficien, the stage hydrograph compare error of stage hydrograph and actual measurement is less under Du family's platform lock of calculating, can think that then this roughness coefficien is institute and asks; Otherwise, should suppose the section roughness coefficien again.

Present embodiment supposition section roughness coefficien is respectively 0.031 and 0.035, as seen from Figure 1, water level is higher under the tame platform lock of Du, roughness coefficien is that the stage hydrograph compare error of stage hydrograph and actual measurement is less under the Du family's platform lock that calculated in 0.031 o'clock, and roughness coefficien is that 0.035 o'clock compare error is bigger; Water level is lower under the tame platform lock of Du, and roughness coefficien is that the stage hydrograph compare error of stage hydrograph and actual measurement is bigger under the Du family's platform lock that calculated in 0.031 o'clock, and roughness coefficien is that 0.035 o'clock compare error is less.Flood road roughness coefficien meets general rule with the characteristics that height reduces that rise of water level under the tame platform lock of Du, can adopt.

Because the calculating section in big vast road has 28, all typings of inconvenience, therefore, only list 4 in Table 1 in and calculated section (promptly being numbered 1,26,27,28), reflection is that the supposition roughness coefficien is the part result of calculation under 0.035 situation in the table one, the time period of extract calculating be October 13 nineteen eighty-three 10 up to 21 o'clock, suppose that roughness coefficien is in the also inconvenient typing table one of result of calculation under 0.031 situation.

Table one: the advance of freshet result of calculation (extraction) during supposition n=0.035

Claims (1)

1. hold the method that increment relation is found the solution holy dimension equation group simulation advance of freshet with groove, it is characterized in that: with St.Venant fundamental differential group after the Δ t period is reduced to water balance equation and kinetic equation formula, can calculate the section groove according to water balance equation (1) and hold increment, according at the beginning of the period, kinetic equation formula (2), the kinetic equation formula (3) at end also can calculate the section groove and hold increment.

Equation group is as follows:

$\mathrm{\&Delta;}{Z}_1=\frac{{{\stackrel{\&OverBar;}{Q}}_1}^2}{{{\stackrel{\&OverBar;}{K}}_1}^2}L---\left(2\right)$

$\mathrm{\&Delta;}{Z}_2=\frac{{{\stackrel{\&OverBar;}{Q}}_2}^2}{{{\stackrel{\&OverBar;}{K}}_2}^2}L---\left(3\right)$

More than in each equation: L is section degree long (m); Δ t is period (s); Z is water level (m); Q is flow (m ³/ s); W is section groove water retention capacity (m ³), Δ W=(W2-W1); Δ Z calculates the water-head (m) of section for the upstream and downstream, section;

Be the section average discharge;

Be section average discharge modulus; Each foot-measuring annotate 1,2 be respectively the period at the beginning of, moment at end.

When known conditions is Q _{Last 1}, Z _{Last 1}, Q _{Following 1}, Z _{Following 1}, Q _{Last 2}And Q _{Following 2}(or Z _{Following 2}) time, can find the solution Z according to the relation that groove holds increment Delta W _{Last 2}And Q _{Following 2}(or Z _{Following 2}):

According to equation (1), as supposition Q _{Following 2}(or Z _{Following 2}) value, just can calculate corresponding section groove and hold increment Delta W (water balance).

According to equation (2), (3), continue to utilize the Q of above-mentioned supposition _{Following 2}(or Z _{Following 2}) value, just can calculate the section groove and hold increment Delta W (power).

Press above-mentioned relation by computer program and regulate Q _{Following 2}(or Z _{Following 2}) value, it is equal substantially that the section groove that dual mode is calculated holds increment Delta W, the Z of this moment _{Last 2}And Q _{Following 2}(or Z _{Following 2}) be solution of equations.

When calculating section number is n, then need to increase n-2 section flow rate calculation relational expression solving equation group.

As the i=1 that is numbered with the upstream initial section, the section numbering toward the downstream then is i=2, i=3 successively ... i=n-2, i=n-1, i=n.

When section numbering 1＜i＜n, according to water balance equation (1), hold increment as the flow of section and section groove and represent with two-dimensional array, then between section i-1 and the section i, the water balance relational expression between section i and the section i+1 is respectively:

Q(i-1，1)+Q(i-1，2)-Q(i，1)-Q(i，2)＝2ΔW(i-1，i)/Δt---------------(4)

Q(i，1)+Q(i，2)-Q(i+1，1)-Q(i+1，2)＝2ΔW(i，i+1)/Δt---------------(5)

Join solve an equation (4), (5):

Q(i，2)＝(ΔW(i，i+1)-ΔW(i-1，i))/Δt+(Q(i-1，2)+Q(i+1，2))/2+(Q(i-1，1)+Q(i+1，1)-2Q(i，1))/2------------------------------(6)

Also can be similar to the accepted way of doing sth forms such as (7) of writing:

Q(i，2)＝(ΔW(i，i+1)-ΔW(i-1，i))/2/Δt+(Q(i-1，2)+Q(i+1，2))/2------(7)

In the formula: first in the two-dimensional array flow Q bracket is the section numbering, second at the beginning of the period Δ t, last code name (just be 1, end be 2); First of holding in the increment Delta W bracket of two-dimensional array groove is upstream, section section numbering, and second be that downstream, section section is numbered.

Therefore, the full section accumulated slot that calculates by the water balance equation holds increment ∑ Δ W (water balance) and is: ∑ Δ W (water balance)=Δ t (Q (1,1)+Q (1,2)-Q (n, 1)-Q (n, 2))/2--------------(8)

The full section accumulated slot that calculates by the power equation holds increment ∑ Δ W (power) and is:

Images