CN101956382B - 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 all is in below the river flood position with main metropolis, 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 confirm 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 types:

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, in the River Basin Hydrology forecasting model, does not adopt as yet.

2, hydrologic method:

Hydrologic method such as famous " Muskingun method " are the fundamental differential group to be reduced to the water balance equation respectively let out equation with being approximately to hold, and 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, hold increment according to also calculating groove with two kinetic equation formulas at end at the beginning of the period.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 through the simulation advance of freshet; The result finds to there is no need to build the Zhi Hongdao project of former report, if implement research report, 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 with hydrologic method with existing hydraulic method, and different places is also arranged.

The place similar with 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 to let out behind the equation to join with the water balance equation separate, and this method is directly water balance equation and kinetic equation formula couplet to be separated.

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 that after the Δ t period was reduced to water balance equation and kinetic equation formula, the relation of holding increment with the section groove joined the group of solving an equation to this method with St.Venant fundamental differential group.

Concrete grammar is following:

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

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

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

More than various in: L is a section length, the m of unit; Δ t is the period, the s of unit; Z is a water level, the m of unit; Q is a flow, the m of unit ³/ s; W is a section groove water retention capacity, the m of unit ³Δ Z is the water-head of upstream and downstream, section calculating section, the m of unit;

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) and formula (2) and formula (3) in order explaining, 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 2ly draw according to water balance equation (1), 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 the process Q that effluents _{Following 1}～Q _{Following 2}The water body that is surrounded.

As shown in Figure 3 according to kinetic equation formula (2) and (3) drafting, 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 known conditions is Q _{Last 1}, Z _{Last 1}, Q _{Following 1}, Z _{Following 1}, Q _{Last 2}And Z _{Following 2}The time, can find the solution Z according to the relation that groove holds increment Delta W _{Last 2}And Z _{Following 2}Or find the solution Z _{Last 2}And Q _{Following 2}:

According to equation (1), like supposition Q _{Following 2}Or Z _{Following 2}Value can calculate corresponding section groove and hold increment Delta W;

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

Press above-mentioned relation through computer program and regulate Q _{Following 2}Or Z _{Following 2}Value, it is equal basically that the section groove that dual mode is calculated holds increment Delta W, the Z of this moment _{Last 2}And Z _{Following 2}Or Z _{Last 2}And Q _{Following 2}Be solution of equations.

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

Like the i=1 that is numbered with upper reaches initial section, the section of downstream numbering 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 like the flow of section and section groove and represent with two-dimensional array, then between section i-1 and the section i and the water balance relational expression between section i and the section i+1 be 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) and (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, and second is first, the last code name of period Δ t, just is 1, and the end is 2; First of holding in the increment Delta W bracket of two-dimensional array groove is the upper reaches, 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

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

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 method that increment relation finds the solution St.Venant equation group simulation advance of freshet according to the figure of kinetic equation formula (2) with (3) drafting for of the present invention with groove.

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 combine Fig. 1 to do further detailed explanation of the present invention, 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 roughness coefficien of the tame platform lock of Du to week during group's section nineteen eighty-three flood.

Shutting out and tame platform lock to week helping the section is that big vast road is about 21km, and the start-stop place is respectively tame Tai Zha of Du and week group, and big vast road is to guide the Han River flood to get into the passage of 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 and gets by flow formula based on the lock upstream and downstream water level of observation.

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 like the roughness coefficien of in table one, supposing;

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

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

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

After accomplishing 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; Can see from Fig. 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, in table one, only listed 4 calculating sections in and promptly be 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 based on water balance equation (1) and hold increment, hold increment based on also calculating the section groove with kinetic equation formula (2) and the kinetic equation formula (3) at end at the beginning of the period;

Equation group is following:

${\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 a section length, the m of unit; Δ t is the period, the s of unit; Z is a water level, the m of unit; Q is a flow, the m of unit ³/ s; W is a section groove water retention capacity, the m of unit ³Δ W=(W2-W1); Δ Z is the water-head of upstream and downstream, section calculating section, the m of unit;

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 known conditions is Q _{Last 1}, Z _{On 1}, Q _{Following 1}, Z _{Following 1}, Q _{Last 2}And Z _{Following 2}The time, can find the solution Z according to the relation that groove holds increment Delta W _{Last 2}And Z _{Following 2}Or find the solution Z _{Last 2}And Q _{Following 2}:

According to equation (1), like supposition Q _{Following 2}Or Z _{Following 2}Value can calculate corresponding section groove and hold increment Delta W;

According to equation (2), (3), continue to utilize the Q of above-mentioned supposition _{Following 2}Or Z _{Following 2}Value also can calculate the section groove and hold increment Delta W;

Press above-mentioned relation through computer program and regulate Q _{Following 2}Or Z _{Following 2}Value, it is equal basically that the section groove that dual mode is calculated holds increment Delta W, the Z of this moment _{Last 2}And Z _{Following 2}Or Z _{Last 2}And Q _{Following 2}Be solution of equations;

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

Like the i=1 that is numbered with upper reaches initial section, the section of downstream numbering 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 like the flow of section and section groove and represent with two-dimensional array, then between section i-1 and the section i and the water balance relational expression between section i and the section i+1 be 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) and (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, and second is first, the last code name of period Δ t, just is 1, and the end is 2; First of holding in the increment Delta W bracket of two-dimensional array groove is the upper reaches, 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 and is:

∑ΔW＝Δ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 and is:

$\mathrm{\&Sigma;\&Delta;W}=\underset{i=2}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;W}(i-1,i)---\left(9\right).$

Images