JPH0786874B2 - Open channel water level / flow rate prediction method - Google Patents

Description

DETAILED DESCRIPTION OF THE INVENTION [Field of Use of the Invention] The present invention relates to an open channel water level / flow rate for predicting the water level and the flow rate at each time along each location along a flow in an open channel (a channel having a free water surface). Regarding prediction method.

[Background of the Invention]

In the case of one-dimensional flow, the descriptive expression of hydraulics in an open channel is as follows when the flow rate and water level are state quantities. (For more information on this item, see, for example, the Hydraulics Official Collection (JSCE).) Where: Q: flow rate (3) H: water level (4) A: cross section of flowing water (5) B: water surface width (6) i: channel bottom slope (7) q: lateral inflow (8) g: gravity Acceleration (9) n: Roughness coefficient (12) x: Downflow distance (place along the flow) (13) t: Time (14).

This equation is a hyperbolic type nonlinear simultaneous first-order partial differential equation with two dependent variables (Q, H) and two independent variables (x, t). According to the characteristic curve theory of partial differential equations, the following simultaneous equations are It can be equivalently converted into an ordinary differential equation.

here, Is.

Equations (15) and (17) are called characteristic curves, and equations (16) and (18) are called characteristic equations.

Conventionally, when solving a characteristic equation on each of these characteristic curves with a digital computer, the distance x and the time t are divided into several grids having a fixed size, and the equations are converted into difference equations at the grid intervals. The fixed-difference lattice type characteristic curve method, which is solved by the above method, has been put to practical use.

Specifically, FIG. 1, as shown in FIG. 2, the upstream side [Delta] x ⁺ j lattice (lattice spacing of the distance direction several respective arbitrary intervals around the grid j, downstream [Delta] x ^- j; j = 1,2, ...
...; Δx ⁺ j = Δx ^- j _-1 ; Δx ^- j = Δx ⁺ j ₊₁ ), and the time direction is divided into evenly spaced grids (grid spacing Δt) at the intersection of the distance grid and the time grid. The state quantity (flow rate and water level) of a certain grid point is finally obtained. The values of the physical quantities of the points (R and S points in the figure) located at the positions narrowed by the distance grid are created by interpolating the values of the grid points at both ends. The lattice spacings Δxj and Δt are the conditions for stability of calculation, here, For all distance grids, a method has been taken so that it is satisfied in advance over the whole time, or based on the calculation results. The flow chart of the calculation step is shown in FIG.

Step 301: Set initial condition and boundary condition of state quantity (flow rate / water level).

Step 302: Spatial (distance) difference grid spacing Δx ^± j, j = 1,2,
... is set.

Step 303: Set the time difference grid interval Δt.

Step 304: Set a calculation point P on the grid j.

Step 305: Calculate the state quantity at the calculation point P on the grid j.

Step 3051: State quantities Q _P , H _P and phenomenon propagation velocity at point P As the first approximation value of the state quantity Q _Q , H _Q and the development propagation velocity at the point Q which is one hour difference on the same grid j To set. At the start of calculation, these values are created from the initial conditions of the state quantity using equations (15), (17), (19), and (20).

Step 3052: Characteristic curve passing through point P2 Is a point that intersects the time difference grid one hour before the time difference grid interval +
-The R point and the S point are set according to the sign.

Step 3053: The state quantities at the points R and S are created by interpolation from the state quantities at the (j-1) point, the point Q, and the (j + 1) point on the same time difference grid.

Step 3054: The state quantity at the point P is calculated by using the equations (15) to (25) which are differentiated. That is, To use.

If these are arranged, the state quantity at the point P is obtained according to the following equation.

(I) When the point P is a point (interior point) other than the boundary (Ii) When point P is on the upstream boundary (a) When the water level is given on the boundary (B) when the flow rate is given on the boundary _{H P = H S + β S} -ξ S (Q P -Q S) (44) (iii) optionally point P is a point on the downstream end boundaries (i) If the water level is given on the boundary, (B) when the flow rate is given on the boundary _{H P = H R + α R} -η R (Q P -Q R) (46) step 306: the value of the state quantity of the point P obtained in the above step 305 It is determined whether or not the difference between these first approximate values is within the allowable range.

Step 307: If the value does not fall within the allowable range, the state quantity of the point P obtained in the above step 305 and the phenomenon propagation velocity obtained from the state quantity are set to the second approximate value of those at the point P, and in the above step 305, The first approximation value is replaced with the second approximation value and steps 305 and 306 are repeated. After that, the approximation value is replaced with a higher order approximation value and the same procedure is repeated.

Step 308: When it is within the allowable range, the state quantity at that time is set to the value at the point P, and the above-mentioned step is performed for the next calculation point.
Steps after 305 are performed. To update the calculation points, first fix the time grid, update the order of spatial grids on that grid, and then update the order of spatial grids on the time grid one hour ahead. To go.

Step 309: After calculating all grid points, output the result and finish.

In the above calculation flow, as an approximate value of the phenomenon propagation velocity at the point P, the average value of the points P and R in the + direction and the average value of the points P and S in the one direction are taken, and it is determined whether or not it is within the allowable range. There is also a method that incorporates some modifications, such as using the phenomenon propagation velocity for judgment, but basically, the spatial difference interval and the time difference interval are fixed in advance, and then the state of each fixed grid point is fixed. It is common in that the characteristic curve defined by the quantity and the phenomenon propagation velocity is determined by iterative convergence calculation.

The conventional fixed difference grid type characteristic curve method described above has the following points.

(1) Difference equation described in discrete form (29) to (46)
Ordinary differential equations described in continuous form based on (15) ~
In order to match (25) with high accuracy, it is necessary to make the spatial difference grid interval Δx ^± j and the time difference grid interval Δt sufficiently small, and the number of grid points increases, which in turn increases the amount of calculation and increases the calculation time. Takes.

The conditions for the stability of the calculation of the difference equation are Equations (26) and (27), and the strictness of the solution as the difference equation is Equations (26) and (2
Since it is the time when the equal sign condition of 7) is satisfied, Δx ^± j and Δt only need to satisfy these conditions, and it seems that there is no need to make them sufficiently small. Does not match unless the following conditional expression, Δx ^± j → 0 (47) Δt → 0 (48), is satisfied at the same time, so the strictness as an ordinary differential equation is always guaranteed even though the strictness as a difference equation can be maintained. It depends on what is not done.

(2) The condition for the stability of the calculation of the difference equation is
(27) is to be satisfied for all time in all spatial grids. Is different for each problem to be solved, and even the same problem changes with time, so Δx ^± j, Δt are determined in advance so that these conditions are satisfied for all problems at all grid points. Is difficult, and it takes a lot of trial and error to make a calculation for each problem several times to decide. Further, depending on the problem, even if this method is used, it will be stable over the entire time Δ
x ^± j, Δt may not be determined.

(3) The strictness of the solution as a difference equation is expressed by equations (26) and (2
When the equality condition of 7) is satisfied, and even if the stability condition is satisfied, if it deviates from the equality condition, the information of points (interpolation points R and S) other than the grid point P is used for calculation. Performance deteriorates, that is, accuracy deteriorates. The degree of deterioration depends on the position of the interpolation point (that is, the phenomenon propagation velocity, Since it depends on and the grid spacing, it differs for each grid point and the uniformity is not guaranteed.

By making Δx ^± j and Δt sufficiently small within the range where stability is satisfied, the uniformity measured by the absolute value can be obtained, but as a result, the number of grid points increases as in the case of (1) above. The amount of calculation becomes enormous and the calculation takes time.

(4) Lattice setting is c ≧ v ≧ 0 (49) When the condition is satisfied, that is, The following conditions are satisfied. This condition is a case where the flow direction is a forward direction and is a normal flow or a limit flow.

Therefore, when the condition that v> c ≧ 0 (52) or v ≦ 0 (53) is satisfied, that is, Or If the above condition is satisfied, the calculation cannot be performed normally as it is. This condition corresponds to the case where the flow direction is the positive flow direction or the flow direction is the reverse direction.

From the above, it cannot be applied to the problems such as the generation of a jet flow, the flow in the opposite direction, and the mixture of a normal flow and a jet flow.

[Object of the Invention]

The object of the present invention is to solve the above-mentioned disadvantages of the conventional method.
By removing the 1st, 2nd, 3rd, and 4th points, the solution of the differential equation is uniformly and highly accurately and stably with a small amount of calculation in comparison with the solution of the original differential equation over all calculation points. In comparison of, uniform and highly accurate and stable with a small amount of calculation,
It is to provide an open channel water level / flow rate prediction method that can calculate even when mixed with superficial flow and reverse flow.

[Outline of Invention]

In order to achieve the above-mentioned object, in the present invention, the difference velocity is adapted every moment during the calculation so that the phenomenon propagation velocity and the difference velocity of the difference equation system are equal, and the adaptation is the spatial difference interval and the time difference. The feature is that both of the intervals are made small at the same time so that the error falls within a predetermined range.

Example of Invention

First, the principle of the present invention will be described below. In the present invention, the physical fact that the stored physical quantity (energy or substance) and the physical phenomenon (change in state quantity such as flow rate or water level) accompanying them are propagated at a finite velocity in the open channel is used. . This finite phenomenon propagation velocity is actually Is.

On the other hand, the stability condition for the calculation of Eqs. (26) and (27) is the phenomenon propagation velocity of the difference equation system. And differential velocity (ratio of spatial grid spacing and temporal grid spacing of difference grid) It shows the relationship with and can be physically interpreted as follows. That is, interpolated points (R point, S point)
When only the information of the grid points is used without using (1), the differential velocity is the velocity at which the phenomenon propagates on the differential lattice when numerically calculating with the differential equation. This means that the structure of the difference equation in which the state quantity at the calculation point is calculated using the state quantity at the adjacent grid point before Δt time is such that the energy or substance at the adjacent grid point is transported to the calculation point at the above-mentioned difference velocity.
It is meant to include the basic assumption of being propagated.

Therefore, (2) If the stability condition is not satisfied, It is calculated that energy and matter are transported and propagated on the difference grid at a velocity smaller than the phenomenon propagation velocity calculated according to the conservation law of the difference equation system (in this case, as shown in Fig. 4, the characteristic curve Either one of the (line with arrow) appears outside the triangle surrounded by the dotted line and the time constant line with the point Q.) As a result, energy and matter are stored in the micro region in the difference equation system. It will break the rules and be accumulated. This is equivalent to artificially giving a high energy density or material density that cannot occur naturally unless a forced external force is applied to the minute region, as a boundary condition for calculating the response amount of the adjacent region. Therefore, in the subsequent calculation of the amount of state of the adjacent region, there is a possibility of calculating a value that cannot occur in an actual phenomenon in order to satisfy the conservation law in the difference equation system. Since this value is used as a boundary condition in the subsequent calculation of the state quantity, it may become a vicious circle, and the calculation may gradually become unstable. At this time, the calculation error is of course accumulated.

(3) If the inequality sign stability condition is satisfied, Energy and matter are transported and propagated (dispersed) on the difference lattice at a speed faster than the phenomenon propagation speed of the difference equation system (a phenomenon that can actually occur if energy and matter density are low) , It shows the effect that the fluctuation becomes smaller than the reality, and the calculation is stable. However, the calculation error is accumulated.

(4) When the stability condition of the equal sign (strictness condition as difference equation) is satisfied, Since energy and matter are transported and propagated on the difference grid at a velocity that matches the phenomenon propagation velocity of the difference equation system, the instability caused by this and the error from the correct solution of the difference equation are eliminated.

Based on the above physical interpretation, it can be said that if the phenomenon propagation velocity and the difference velocity of the difference equation system are equal, the stability and rigor of the calculation as the difference equation are guaranteed.

The method of the present invention realizes an adaptive differential (lattice) type characteristic curve method in which the differential velocity is adapted every moment during the calculation so that the phenomenon propagation velocity of the differential equation system becomes equal to the differential velocity. This is to make it possible to obtain a solution that has stability and rigor in the calculation of the difference equation. Differential speed is Therefore, the method of adaptation is as follows: spatial difference adaptive type that changes Δx ^± j, time difference adaptive type that changes Δt (above, see Japanese Patent Application No. 58-40679), and Δx ^± j and Δt. There is a space-time difference adaptive type that changes both, and the present invention relates to the last space-time difference adaptive type.

The principle of the present invention will be described in more detail.

To match the differential equation and the difference equation, formally, while satisfying the stability strictness constraint condition of the difference equation with respect to the spatial difference interval Δx ^± j and the temporal difference interval Δt,
It has already been stated that Δx ^± j and Δt should be close to 0 at the same time. This is because the order of approximation error is expressed by the function of Δx ^± j and Δt, regardless of the approximation method of the difference that transforms the equation from continuous form to discrete form, such as first-order, second-order, etc. However, when Δx ^± j → 0 and Δt → 0, the function value becomes 0. By the way, in the numerical calculation by the digital computer, Δx ^± j and Δt cannot be set to 0, and the difference must be finite. Then, there arises a problem of how much Δx ^± j and Δt can be reduced to obtain sufficient matching. Japanese Patent Application Sho 58-40
In No. 679, it was stated that the adaptive difference method satisfying the requirements of calculation stability and strictness as a difference equation can generate only the minimum necessary calculation points in time if Δx ^± j is given in advance. . This means that when Δx ^± j is determined, Δt is inevitably determined at each time and each place (that is, the resolution with respect to time is also determined). Means that there is one degree of freedom left regarding how to determine Δx ^± j. We believe that this degree of freedom has the key to solving the problem of conformity.

Considering the agreement between the differential equation and the difference equation from the viewpoint of the behavior of the solution in the space-time domain, as shown in FIG. 5, a physical phenomenon described in a continuous form in that domain (actual physical phenomenon Note that the curved surface (the curved surface shown by the solid line) 50 that represents the state quantity that expresses There is a curved surface 53 (curved surface indicated by a broken line) formed by interpolating (interpolating / extrapolating) the value of the state quantity at (●) at a distance 54 within the allowable error. Is.

In FIG. 5, an arrow 55 on each curved surface represents a state quantity change vector.

Here, if the interpolation formula is well formed (that is, if it is faithful to the described physical phenomenon), the density of calculation points may be low, but in practice such an interpolation formula cannot be created, Instead of lowering the interpolation formula, the density of calculation points is increased. On the other hand, if only one kind of interpolation formula is used in the region, the density of the calculation points is high at the space-time point (601 in the figure) where the state quantity change curvature is large, and the change curvature is small, as shown in FIG. However (at 602 in the figure), it is desirable to lower the value because it reduces variations in interpolation error.

The adaptive difference method is not a change of the state quantity itself, but the density of the calculation points is low at the space-time point where the phenomenon propagation velocity is large and the density of the calculation point is high at the small space-time point. It is adapted to automatically adjust the axial density. However, in the adaptive difference method Δ
Since x ^± j is fixed without depending on the phenomenon, the variation of the error as the difference equation is reduced, but the consistency with the original differential equation causes an error in the interpolation in the spatial direction. Depends on the size of Δx ^± j. On the contrary, if the allowable error is specified, it means that the calculation points that satisfy the allowable error, that is, the methods of determining Δx ^± j and Δt can be created in consideration of the interpolation method.
Further, it means that if the allowable error is determined, the temporal and spatial resolution is inevitably determined.

In the characteristic curve method, it is considered that the phenomenon propagates along the characteristic curve and the state quantity changes. Therefore, in order for the state quantity surface of the continuous form system to agree with the state quantity surface of the discrete form system, 2
The characteristic curves of the two systems must match.

In the differential characteristic curve method, as shown in FIG. 7, the characteristic curve (701,703) passing through each calculation point is regarded as a straight line (702,704), and the point (705,706) necessary for calculation is calculated on the straight line by the convergence calculation. Ask. This is equivalent to piecewise linear approximation of the characteristic curve. At that time, the value of the state quantity of those necessary points is created by linear interpolation from the value of the state quantity of the adjacent calculation points (707,708,709), and the state quantity of the calculation point is calculated using the created value. Therefore, in this process, when performing piecewise linear approximation of the same characteristic curve, and between characteristic curves (between adjacent calculation points)
An error occurs when the value is interpolated by. The former can be regarded mainly as an error directly related to the consistency between the continuous form and the discrete form, and the latter can be regarded as an error mainly related to the consistency indirectly through the uniformity of the interpolation error between the characteristic curves.
If these errors can be reduced, it can be said that the discrete solution can be made closer to the continuous solution. The present invention relates to a method of reducing these errors. In order to reduce the error caused by the piecewise linear approximation of the same characteristic curve, generally, when the curve is piecewise linearly approximated, the approximation error can be reduced by dividing the curve as finely as possible. At that time, if the distance is divided into equal intervals, the approximation error differs depending on the place, and in order to make it less than the allowable error, there are places that need to be divided unnecessarily finely, which leads to an increase in the amount of calculation. Therefore, in order to equalize the error, a method in which a gentle change in the slope of the curve is coarse and a sharp change is finely divided can be considered. The present invention also employs this idea, but the challenge is to do this for a continuous type characteristic curve of unknown shape.

Now, as a precondition, if Δx ^± j is made smaller, the adaptive difference method is used so that the characteristic curve in the discrete form matches that in the continuous form. Moreover, as shown in FIG. 8, if the characteristic curve can be linearly approximated within the allowable error, Δx
^{If ±} j (Δx 'in the figure) is given, smaller Δx ^± j (Δx in the figure) can also be linearly approximated within the tolerance,
And the slope of the characteristic curve in both cases (801, 802 in the figure) there is little difference, and they differ only in the range of the tolerance with the continuous type gradient (803 in the figure). From these two things, if Δx ^± j is gradually increased from a small value that can be said to be consistent with the continuous form, It can be said that the consistency with the continuous form is maintained within the error range as long as it does not deviate from the value when the value is small by more than the allowable error, and even the characteristic curve of the unknown continuous form can be piecewise linearly approximated within the allowable error. become.

On the other hand, if the continuous type and the discrete type characteristic curves match within the allowable error range, Is small for the same Δx ^± j at the point where the slope of the characteristic curve is gradual, and conversely is large at the point where the slope is abrupt, so that between two adjacent points on the same characteristic curve in the discrete form If Δx ^± j is adjusted to various values so that the difference in Δ is within the permissible error at any point, the minimum required Δx ^± j according to the steepness of the gradient change can be found.

Therefore, if the initial value of Δx ^± j is set to a sufficiently small value in the numerical calculation and this adjustment is constantly performed during the subsequent calculations, the continuous and discrete characteristic curves must have the minimum required Δ
A piecewise linear approximation can be made with x ^± j and Δt.

The absorption calculation is used to adjust Δx ^± j.
If the adjusted value of the neighboring points is used as the initial value, as long as the slope of the characteristic curve does not change rapidly. As long as the rate of change of does not change significantly), there is a high probability that convergence will be faster, and an increase in the amount of calculation can be suppressed.

Taking the above into consideration, specifically, the error caused by the piecewise linear approximation is reduced by the following procedure.

In the following (A1) to (A7), the calculation point P defined by certain Δ ▲ x ⁺ _P ▼ and various quantities related to the point P are calculated by convergence calculation. After (A8), Δ ▲ x ⁺ _P ▼ Regarding finding the minimum required size.

In FIG. 9, it is assumed that the state quantity at (A1) t = t ₀ is given or can be calculated by interpolation.

(A2) The phenomenon propagation velocity at t = t ₀ can be calculated from the state quantity.

(A3) The minimum required Δx ^± j at t = t ₀ is given or can be calculated by interpolation with an accuracy within the allowable error.

(A4) Calculation point P on the characteristic curve that passes through point R (x ₀ , t ₀ ).
To determine the ^{_{(x 0 + Δ ▲ x +}} P ▼, t 0 + Δt) Δ ▲ x + P
As the initial value of ▼, take Δx ^± j at the proximity calculation point.

(A5) Phenomenon propagation speed at point R A straight line that has a gradient of and passes through the point R and a straight line x = x ₀ + Δ ▲ x ⁺ _P
Let the intersection with ▼ be the first approximation of the point P, and set the point Q (x ₀ + Δ ▲ x ⁺ _P
, T ₀ ) as the first approximation of the state quantity at point P
Phenomenon propagation speed Compute a first approximation of

(A6) Point P The first approximation of the point R of (A5) Instead of, make a second approximation of point P on x = x ₀ + Δ ▲ x ⁺ _P ▼, and The first approximation of the point S is made on t = t ₀ from the first approximation of
A second approximation of the state quantity at the point P is calculated from the state quantity at the point S, and the phenomenon propagation velocity is calculated using the second approximation. The second approximation of is calculated.

(A7) A similar method is used to improve the approximation and converge the value related to point P, and the state quantity and phenomenon propagation speed at that time And differential speed Be the solution in the discrete form with respect to the initial value of the spatial difference interval Δ ▲ x ⁺ _P ▼.

(A8) Next, as shown in FIG. 10, Δ ▲ x ⁺ _P ▼ the initial value per case of less than (corresponding to Δ ▲ x ⁺ _P ▼ Figure midpoint P _9), similar to the above Method to find a discrete solution (corresponding to point P ₃ in the figure). How small ΔΔx ⁺ _P ▼ is to be set separately in consideration of convergence acceleration.

(A9) State quantity obtained in (A7) and (A8) above and /
Alternatively, it is determined whether or not the difference in phenomenon propagation speed is larger than the allowable error range.

(A10) If it is larger than the allowable error range in (A9) above,
It is considered that the two points (points P ₉ and P _{3 in the} figure) of the point R and the point P corresponding to the small Δ ▲ x ⁺ _P ▼ do not exist on the same characteristic curve, and ΔΔ is further increased than in the case of (A8) above. x ⁺ _P ▼ is reduced (Δ ▲ x ⁺ _P ▼ corresponding to the point P ₂ in the figure), and this is repeated until it is within the allowable error (points P ₃ and P _{2 in the} figure).

(A11) If it becomes smaller than the allowable error range in (A9) above, two points, point R and point P (point P _{2 in the} figure) corresponding to a small Δ ▲ x ⁺ _P ▼, are on the same characteristic curve. , And the continuous solution is regarded as approximating within the allowable error range, and Δ
▲ x ⁺ _P ▼ is larger than these two cases (Δ ▲ x ⁺ _P ▼ corresponding to points P ₂ and P ₃ in the figure) (Δ ▲ corresponding to point P ₄ in the figure).
x ⁺ _P ▼), a discrete form solution is obtained by the same method as described above, and the same comparison and determination as in (A9) is performed. This is repeated until the difference exceeding the allowable error range (corresponding to point P _{5 in the} figure) appears (corresponding to point P _{4 in the} figure), and ΔΔx ⁺ _P ▼ is continuously increased. Δ ▲
How much x ⁺ _P ▼ is increased at one time is determined separately in consideration of the convergence acceleration.

(A12) From the above, the minimum Δ ▲ x ⁺ _P ▼ (corresponding to point P _{4 in the} figure) and Δ ▲ x ^- _P ▼ required to obtain a discrete-form solution that matches the continuous form within the allowable error range are determined. .

(A13) In (A11) above, three or more Δ ▲ x ⁺ _P ▼
(Equivalent to points P ₄ , P ₃ , P _{2 in the} figure) is compared and discriminated, and when three or more points are continuously within the allowable error, the minimum Δ ▲ x ⁺ _P ▼
(Figure midpoint P ₂ or equivalent) corresponding points P (Figure midpoint P ₂₎ and the point R are on the same characteristic curve actually adopted is a method for strict conditions so on. This reduces the risk of adopting a solution that satisfies the discrete form (difference equation) but not the continuous form (differential equation).

(A14) In the above (A4), at the start of all calculations, a small value is provided as the initial value of Δ ▲ x ⁺ _P ▼ such that the continuous form and the discrete form are matched.

In (A9) of the above procedure, not only the comparison with the permissible error is performed by the phenomenon propagation velocity indicating the slope of the characteristic curve,
Although it is possible to use state quantities, this is because the final evaluation of the error in numerical calculation is done with respect to the dependent variable or state quantity in the original differential equation, as described above. This is to think that it should be done.

In FIG. 9, all black dots 91 indicate calculation points.

Reducing interpolation error between characteristic curves Generally, when linearly interpolating the value of a point between curves from the value of a point on the curve, the density of the curve should be as high as possible to reduce the interpolation error. Good. At that time, as shown in FIG. 11, the value (4 ′) of an arbitrary point (point 1104) inside a triangle composed of three adjacent points (points 1101, 1102, 1103) on at least two adjacent curves is defined as , The difference (1105) between the value and the value (4 ″) calculated by linear interpolation from the values (1 ′, 2 ′, 3 ′) at the three adjacent points is within the allowable error range. Any density will do.

Also in the case of the characteristic curve, a method of generating the adjacent characteristic curve (adjacent calculation points) that satisfies the requirement given by the above-mentioned tolerance will be described below with reference to FIG.

Now, referring to FIG. 12, as a precondition, the same characteristic curve is piecewise linearly approximated within the allowable error range by using the method described above. It is assumed that the state quantity of the inner point of the triangle PRS surrounded by and t = t ₀ can be calculated within the allowable error range by linear interpolation. The problem here is that when a new calculation point P'is created outside the triangle PRS, an arbitrary point inside the triangle consisting of the point P'and at least one point in the vicinity of the point P'is known to be on another characteristic curve. The point values (phenomenon propagation velocity and / or state quantity) are to be matched with the interpolated values from the values of the three endpoints of the triangle within the allowable error range. When making a new calculation point P ′ outside the triangle PRS, if it is desired to exist on the same characteristic curve as the point P, the point R ′ relating to the point P ′ should be placed on the line segment RP. As mentioned before. Here, a point R'is placed on a point other than the line segment PR in the triangle PRS, and the point P'is changed to the point P.
Consider a case where it is placed on a characteristic curve different from.

Specifically, the following procedure reduces the interpolation error between characteristic curves.

(B1) As the time line (t = t ₁ ) at which point R ′ exists, select the time line at which the most delayed pre-calculated point exists at that time, and the state quantity is given on that line. .

(B2) From the point R ′ near the line segment PS in the triangle PRS to the point P
At a position with the spatial difference interval Δ ▲ x ⁺ _P ▼ of which is the initial value, a point P ′ existing on the same characteristic curve as the point R ′ and a point S ′ on the other characteristic curve passing through the point P ′. Obtained according to the method described above. It is considered that the value of an arbitrary point in the newly formed triangle P'R'S 'can be calculated from the values of the three end points by interpolation within the allowable error range.

(B3) An arbitrary point (for example, the point P on the line segment P'P) is set inside the triangle surrounded by the line segments P'P, PS, P'R ', and the value of the point is the point R'. Spatial adaptive difference method for determining Δx ^± j by fixing Δt from the value on t = t ₁ (Japanese Patent Application No. 58-4067).
(See No. 9).

(B4) It is checked whether or not the calculated value matches the value created by linear interpolation from the values of the end points of the triangle within an allowable error range.

(B5) If they match, the point P'is regarded as a point that satisfies the requirement of the allowable error.

(B6) If they do not match, t = t ₁ upper point R at a point in the triangle PRS
The point R'is placed on the opposite side of the point R ', and the point P'
The point P'is created by the same method as that used for creating.

(B7) The same test as that performed for the point P ′ and the point P ′ is performed, and these steps are repeated until a match is found within the allowable error range.

(B8) If the points P'or P ', etc. are found to be in agreement, the processing will be described in (B9) by taking the case of the point P'as an example.

(B9) point P 'on the opposite side to the time-space and the point P to the boundary of is in the t> t _1, and triangular P'r'S' nearest the P 'point on the outside of When the calculated point is P ″, an arbitrary point (for example, a point on the line segment P′P ″) is set inside the triangle surrounded by the line segments P′P ″, P ′S ′ and P ″ R ″. At that point, the same test as that performed at the point P is performed, and a point similar to the point P'is created between the points P'and P "until the condition of the tolerance is satisfied.

After the above steps are completed, regarding the points P, P ″, and the newly created points P ′ and P ′, if any two adjacent points are taken out and considered, on the time line passing through the points that are delayed in time. The value of can be calculated by linear interpolation within the allowable error range between the intersection of the characteristic curve passing through the point advanced in time and its time line and the point delayed in time.

(B10) Return to step (B1) and proceed to the next calculation point.

(B11) If t ₁ = 0, first create a calculation point P,
Next, a point P ′ may be created and the same procedure may be performed. However, since there is no point corresponding to the point P ″, a point corresponding to the point P ″ is created in the same manner as the point P ′.

Next, an embodiment of the present invention will be described with reference to FIG.

An open channel water level / flow rate prediction method according to the present invention includes a state quantity (flow rate / water level) initial condition / boundary condition setting unit 1, a spatial difference interval initial trial value setting unit 2, a calculation reference origin R setting / updating unit 3, Calculation point P'creation part 4 with small interpolation error between characteristic curves, calculation end judgment part 5, output part 6, calculation point P'creation part 41 with small error due to piecewise linear approximation of characteristic curve, convergence judgment part A 42, a changing unit 43 of the calculation reference origin R ′, a creating unit 411 of a point P _i that satisfies the stability strictness of the difference form, a convergence determining unit B 412, a spatial difference interval changing unit 413, a state quantity of the point P _i , etc. The calculation unit 4111, the time difference interval adaptation unit 4112, and the convergence determination unit C 4113 are included.

The operation of the method according to the invention is described below.

First, (C1) state quantity (flow rate / water level) initial condition / boundary condition setting unit 1 sets initial conditions and boundary conditions in the time / space domain to be calculated. In (C2) Spatial difference interval initial trial value setting unit 2, the spatial difference interval Δx ⁺
Set the initial trial value for _j .

Then, using the initial trial value of (C3) Δx ⁺ _j , at 41,
A point P on the characteristic curve passing through the origin (point R) and a point S on another characteristic curve passing through the point P, which satisfy the stability strictness of the differential form within the range of the allowable error, are obtained. This is a Δx ⁺
A solution satisfying the stable strictness of the difference form with respect to _j is obtained by a convergence calculation by the method of adapting the time difference interval at 411, and the decision result at 412 is checked so that the solution converges to the solution of the differential form. In the method, the spatial difference interval Δx ⁺ j is changed and the calculation is repeated.

In the method of adapting the time difference interval in the above 411, the state quantity and the phenomenon propagation velocity are calculated at 4111 for a given point P _i , and the time difference interval Δt is adapted at 4112 to the phenomenon propagation velocity, and the result is obtained. , 4113 determines that the state quantity and the phenomenon propagation speed have converged, and then passes control.

By (C3) above, (C4) the spatial difference interval Δx ^± j and the time difference interval Δt are adaptively adjusted so as to satisfy the stability strictness of the differential form within the allowable error range.

Next, point (C5) S is set as point R for calculating a new calculation point, and Δx ⁺ j adjusted in (C4) above is set as an initial value (C
Similar to 3), find points P and S of the new calculation points.

(C6) Repeat the above (C5) until the end point in the spatial direction comes.

Next, (C7) Create new calculation points between each calculation point obtained up to (C6) so as to satisfy the requirement from the allowable error regarding the interpolation error between the characteristic curves in 41 by the above method. To do. This is because the method of (C3) is repeated while the convergence reference is made at 42 with respect to the calculation reference origin R'changed at 43.

Next, (C8) a point which is the latest in time and is not covered by the characteristic curve from the adjacent calculation point is searched for, and the point is assumed to be on the time line where the point R ′ exists, and the method described above is used. Create a new calculation point in 3.

At this time, in some cases, the characteristic curve of the new calculation point may be made to cover the adjacent calculation points and the calculation points outside thereof, which has the effect of reducing the number of calculation points.

Next, (C9) The above (C8) is repeated until the calculation of all the calculation areas is completed, the judgment is made at 5, and the result is outputted at 6.

〔The invention's effect〕

According to the present invention, (1) the time and space difference intervals are set at each place and each time so that the condition of the stable exact solution of the difference equation is satisfied within the range of the interpolation error of the points R and S. To adapt to
Stable, uniform and highly accurate calculation can be realized over all calculation points.

(2) Phenomenon propagation speed Since there is no restriction on the magnitude relationship and the positive / negative relationship of, the normal flow, the limit flow, the forward and reverse directions of the superficial flow and the flow, and their mixed state can be accurately calculated.

(3) Fixed-difference grid method in which the entire calculation area is covered with a fine grid in order to satisfy the stability condition in order to apply the minimum necessary time and space difference according to the change in the phenomenon (the magnitude of the phenomenon propagation speed). Compared with, the number of calculation points is significantly reduced, and it can be applied to online calculation control by a microcomputer.

[Brief description of drawings]

FIG. 1 is a schematic explanatory view of a fixed difference grid type characteristic curve method, and FIG.
The figure is a detailed explanatory diagram of the fixed difference grid type characteristic curve method. Fig. 3 is a flow chart of calculation of the fixed difference grid type characteristic curve method. Fig. 4 is an example when the calculation becomes unstable. Fig. 5 is a continuous type. And FIG. 6 are explanatory diagrams of matching between the discrete type and the discrete type, FIG. 6 is an explanatory diagram of equalizing the interpolation error, FIG. 7 is an explanatory diagram of the error in the differential type characteristic curve method, and FIG. 8 is a continuous type characteristic curve and the discrete type. FIG. 9 is an explanatory diagram of coincidence of characteristic curves, FIG. 9 is an explanatory diagram of calculation corresponding to a given spatial difference interval, FIG. 10 is an explanatory diagram of calculation of a minimum required spatial difference interval, and FIG. 11 is an arbitrary curve interval. FIG. 12 is an explanatory diagram of interpolation of a value on a dot, FIG. 12 is an explanatory diagram of a method of reducing an interpolation error between characteristic curves, and FIG. 13 is an embodiment of the present invention.

Claims (1)

[Claims] 1. A step of measuring the water level or flow rate of an open channel, a step of setting the water level or flow rate obtained in the above measuring step as boundary conditions and initial conditions, and a difference equation for describing the state of the water level and flow rate of the open channel. Between the step of setting the initial trial value of the spatial difference grid interval, the step of setting / updating the origin of time and the reference point for the prediction calculation of the water level and discharge state, and the characteristic curve used for the prediction calculation of the water level and discharge The first creation step of creating a calculation point of a place and time where the interpolation error is small, the step of determining the end of the prediction calculation, the step of outputting the predicted value of the water level and the flow rate of the result of the determination
And a second creating step for creating a calculation point at a place and time where a calculation error caused by a piecewise linear approximation of a characteristic curve used for predictive calculation of water level and discharge is small. The first convergence determination step for determining whether or not the calculation error due to the linear approximation has converged to a small value, and when the convergence does not occur in the first convergence determination step, the place to be the reference for predictive calculation and the origin of the time are changed. The second creating step is a third creating a point at a place and time that satisfies the stability strictness of the difference equation, which is a discrete physical model that describes the state of water level and discharge. Step, a second convergence determination step of determining whether the above-mentioned stability rigor is satisfied, and a spatial difference grid interval of the above difference equation is changed if the convergence does not occur in the second convergence determination step. The step of, consist city, the third creation step, the water level of calculation points to predict calculate water levels and flow rates, flow rates, transport and
Calculating the propagation velocity, from the relationship between the calculated transport / propagation velocity and the spatial difference grid interval, the ratio of the spatial difference grid interval and the temporal difference interval becomes the absolute value of the transport / propagation speed. The steps of adaptively calculating the temporal difference grid interval so that they are equal to or larger than that, determining whether or not the water level and flow rate and transport / propagation velocity have converged after passing through the above steps, and converging If not, a third convergence determination step of passing control to the step of recalculating the water level, the flow rate, and the transportation / propagation velocity using the calculated temporal difference grid, and an open channel. Water level / flow rate prediction method.