CN105046069A - Method for predicting numerical value of water temperature in reservoir of large hydropower station - Google Patents

Abstract

The invention relates to the field of water conservancy and hydropower engineering, and in particular relates to a method for predicting a numerical value of water temperature in a reservoir of a large hydropower station. The method comprises the following steps: (1), collecting basic parameters of a hydropower station damsite; (2), establishing a numerical analysis model according to the collected basis parameters; (3), determining the collected basis parameters as low-flow year parameters, median-flow year parameters and high-flow year parameters according to a low-flow year, a median-water year and a high-flow year; (4), numerically simulating a water temperature change rule under laminar flow heat exchange of the reservoir within a year respectively according to the low-flow year parameters, the median-flow year parameters and the high-flow year parameters; (5), obtaining a low-value result, a median-value result and a high-value result of water temperature of the reservoir in the front of the dam at the upstream during the operation period of the reservoir; and (6), determining a distribution rule and a quantity value of water temperature of the reservoir during the operation period so as to obtain an envelope type prediction result of the water temperature of the reservoir. According to the invention, influence factors, such as the shape of the reservoir, the hydrology-weather condition in a reservoir area, the operation condition of the reservoir and the initial water storage condition of the reservoir, are considered sufficiently, so that the predication precision is greatly improved.

Description

A kind of large hydropower station storehouse coolant-temperature gage Numerical Predicting Method

Technical field

The present invention relates to Hydraulic and Hydro-Power Engineering field, particularly relate to a kind of large hydropower station storehouse coolant-temperature gage Numerical Predicting Method.

Background technology

Water temperature of reservoir is the important temperature boundary condition of of power station concrete dam, is one of dam temperature stress and temperature controlled important factor in order.Upper pond water temperature will directly have influence on the distribution in dam runtime stable (quasi-steady) temperature field, particularly for dam body foundation confining region, bottom of the reservior, upstream water temperature will directly have influence on the basic temperature difference and the temperature control criterion of dam, thus have influence on the Design of Temperature Control of concrete dam.How before reservoir builds up, reasonably prediction reservoir build up after distribution of water temperature situation, be an important prerequisite of Dam Designs in Last.

The reservoir of large hydropower station is by after retaining Cheng Ku, the change of water temperature is a very complicated phenomenon, by the control of many factors, major influence factors has four aspects: the shape of reservoir, the initial water storage condition of reservoir area hydrometeorology condition, Reservoir Operation Conditions and reservoir.In existing correlation technique, build up the Forecasting Methodology of Ba Qian storehouse, rear upstream coolant-temperature gage for reservoir, except partly can considering the influence factor of weather and regional condition, the initial retaining situation of the shape of concrete reservoir, Reservoir Operation Conditions and reservoir cannot be considered.In hydrometeorology condition, also cannot consider the factors of influence such as wind speed, cloud amount, water surface evaporation, river water temperature, reservoir inflow, warehouse-in silt content.Therefore for the different reservoirs of identical region, the water temperature of reservoir utilizing the method for prior art to predict out is substantially identical, cannot reflect the characteristic of concrete reservoir, and the accuracy therefore predicted is not high.

Summary of the invention

The object of the present invention is to provide a kind of large hydropower station storehouse coolant-temperature gage Numerical Predicting Method, comprise the steps:

Step 1, gathers the underlying parameter of Dam Site; Wherein, underlying parameter comprises the initial retaining parameter of reservoir form parameter, hydrometeorological parameter, reservoir operational factor and reservoir;

Step 2, the underlying parameter according to gathering sets up numerical analysis model;

Step 3, is defined as low flow year parameter, normal flow year parameter and high flow year data according to low flow year, normal flow year and high flow year by the underlying parameter of collection;

Step 4, according to low flow year parameter, normal flow year parameter and high flow year parameter, the water temperature Changing Pattern respectively in numerical simulation year under reservoir laminar flow exchange heat;

Step 5, again respectively to not considering upstream Ba Qian bottom of the reservior slag and alluvial, the alluvial of Ba Qian bottom of the reservior, consideration upstream, considering three kinds of situations such as upstream Ba Qian bottom of the reservior slag and alluvial, carry out numerical simulation, obtain low value result, median result, the high level result of reservoir runtime upstream Ba Qian storehouse coolant-temperature gage;

Step 6, according to the low value result of reservoir runtime upstream Ba Qian storehouse coolant-temperature gage obtained, median result and high level result, determines the regularity of distribution and the value of runtime water temperature of reservoir, obtains water temperature of reservoir envelop-type and predicts the outcome.

Further, reservoir form parameter comprises relation between reservoir capacity, the reservoir degree of depth, reservoir level ~ storage capacity ~ storehouse length ~ area; Described hydrometeorological parameter comprises temperature, solar radiation, wind speed, cloud amount, evaporation capacity, reservoir inflow, warehouse-in water temperature, river load content, warehouse-in suspended load; Reservoir operational factor comprises dam safety evaluation mode, irrigating gate position, power station and diversion ability, flash buildings position and discharging capacity, the traffic control situation of reservoir, storage outflow, reservoir level change; Slag situation, upstream cofferdam disposition before ground temperature, first filling temperature, reservoir filling speed, dam when the initial retaining parameter of reservoir comprises first filling season, first filling.

Further, described step 2 specifically comprises: set up model of reservoir according to reservoir form parameter; Set up water body heat conduction model; Set up the heat transfer model of warehouse-in, outbound water body.

Compared with prior art the invention has the beneficial effects as follows: for different reservoir feature separately, the distribution of water temperature rule that envelop-type forecast analysis Reservoir is later, take into full account the influence factor such as the shape of reservoir, the initial retaining situation of reservoir area hydrometeorology condition, Reservoir Operation Conditions and reservoir, substantially increase the accuracy of prediction.

Accompanying drawing explanation

Fig. 1 is the process flow diagram of a kind of large hydropower station storehouse of the present invention coolant-temperature gage Numerical Predicting Method;

Fig. 2 is JP-1 reservoir level ~ storage-capacity curve figure;

Fig. 3 is JP-1 reservoir level ~ curve of areas figure;

Fig. 4 is the long curve map in JP-1 reservoir level ~ storehouse;

Fig. 5 is water temperature of reservoir computation model sketch;

Fig. 6 is JP-1 reservoir annual distribution of water temperature figure (not considering front bottom of the reservior, dam slag);

Fig. 7 is JP-1 reservoir monthly average distribution of water temperature figure (not considering front bottom of the reservior, dam slag);

Fig. 8 is JP-1 reservoir annual distribution of water temperature figure (considering front bottom of the reservior, dam slag);

Fig. 9 is JP-1 reservoir monthly average distribution of water temperature figure (considering front bottom of the reservior, dam slag).

Embodiment

Also by reference to the accompanying drawings the present invention is described in further detail below by specific embodiment.

Shown in ginseng Fig. 1, Fig. 1 is the process flow diagram of a kind of large hydropower station storehouse of the present invention coolant-temperature gage Numerical Predicting Method.

A kind of large hydropower station storehouse coolant-temperature gage Numerical Predicting Method that the present embodiment provides, comprises the steps:

Step S101, gathers the underlying parameter of Dam Site; This underlying parameter comprises the initial retaining parameter of reservoir form parameter, hydrometeorological parameter, reservoir operational factor and reservoir;

Step S102, the underlying parameter according to gathering sets up numerical analysis model;

Step S103, is defined as low flow year parameter, normal flow year parameter and high flow year parameter according to low flow year, normal flow year and high flow year by the underlying parameter of collection;

Step S104, according to low flow year parameter, normal flow year parameter and high flow year parameter, the water temperature Changing Pattern respectively in numerical simulation year under reservoir laminar flow exchange heat;

Step S105, again respectively for not considering upstream Ba Qian bottom of the reservior slag and alluvial, the alluvial of Ba Qian bottom of the reservior, consideration upstream, considering three kinds of situations such as upstream Ba Qian bottom of the reservior slag and alluvial, carry out numerical simulation, obtain low value result, median result, the high level result of reservoir runtime upstream Ba Qian storehouse coolant-temperature gage;

Step S106, according to the low value result of reservoir runtime upstream Ba Qian storehouse coolant-temperature gage obtained, median result and high level result, determines the regularity of distribution and the value of runtime water temperature of reservoir, obtains water temperature of reservoir envelop-type and predicts the outcome.

The Numerical Predicting Method that the present embodiment provides, for different reservoir feature separately, the distribution of water temperature rule that forecast analysis Reservoir is later, take into full account the influence factor such as the shape of reservoir, the initial retaining situation of reservoir area hydrometeorology condition, Reservoir Operation Conditions and reservoir, substantially increase the accuracy of prediction.

In the present embodiment, reservoir form parameter comprises relation between reservoir capacity, the reservoir degree of depth, reservoir level ~ storage capacity ~ storehouse length ~ area; Hydrometeorology parameter comprises temperature, solar radiation, wind speed, cloud amount, evaporation capacity, reservoir inflow, warehouse-in water temperature, river load content, warehouse-in suspended load; Reservoir operational factor comprises dam safety evaluation mode, irrigating gate position, power station and diversion ability, flash buildings position and discharging capacity, the traffic control situation of reservoir, storage outflow, reservoir level change; Slag situation, upstream cofferdam disposition before ground temperature, first filling temperature, reservoir filling speed, dam when the initial retaining parameter of reservoir comprises first filling season, first filling.

In the present embodiment, initial modeling data comprises:

(1) reservoir level ~ storage capacity ~ storehouse length ~ curve of areas;

(2) power station irrigating gate number, size, center line elevation and maximum diversion ability (rice ³/ second);

(3) flash buildings number, size, center line elevation and maximum discharging capacity (rice ³/ second);

(4) for many years on average (or Typical Year, lower same) by ten days, moon temperature (DEG C);

(5) average by ten days, moon river water temperature (DEG C) for many years;

(6) average wind speed (meter per second) month by month for many years;

(7) for many years average solar radiation quantity month by month (kilocalorie/centimetre ². sky);

(8) average cloud amount (one-tenth number) month by month for many years;

(9) average evaporation capacity (evaporation from water surface field, kilograms per centimeter month by month for many years ². sky);

(10) on average suspended load or silt content (kg/m is put in storage month by month for many years ³);

(11) average reservoir inflow (rice month by month for many years ³/ second);

(12) storage outflow (rice is month by month designed ³/ second) and Reservoir Operation Conditions month by month;

(13) upstream reservoir level and downstream tailwater level month by month in design and operation year;

(14) slag situation, upstream cofferdam disposition etc. before ground temperature, first filling temperature, reservoir filling speed, dam when first filling season, first filling;

(15) before reservoir filling, slag elevation before dam, upstream, slag season, slag body material property, slag technique etc.;

(16) dam safety evaluation performance;

(17) water diffusion coefficient; Water Convention diffusion coefficient; The absorption coefficient of water.

In the present embodiment, the establishment step of numerical analysis model comprises:

1) model of reservoir is set up according to reservoir form parameter:

For JP-1 reservoir.Its reservoir form parameter is as follows:

Table 1JP-1 reservoir level ~ storage capacity ~ area ~ storehouse is long

Reservoir level (m)	Volume (hundred million m 3)	Area (km 2)	Storehouse long (km)
				1580	0.000	0.000	3.800
1600	0.005	0.352	12.500
				1640	1.000	2.200	23.300
1650	1.400	3.100	25.100
				1660	2.100	3.724	27.100
1670	2.500	5.100	29.000
				1680	3.300	6.300	31.000
1690	4.100	8.100	33.100
				1700	5.400	10.300	34.700
1710	7.000	12.100	36.200
				1720	8.300	14.100	37.800
1730	10.100	17.200	39.200
				1740	11.800	20.300	41.000
1750	13.800	24.520	42.400
				1760	16.100	27.740	43.500
1770	18.800	31.320	44.600
				1780	21.300	35.030	45.600
1790	24.600	38.690	46.700
				1800	28.540	42.740	47.900
1810	33.010	46.630	49.100
				1820	37.880	50.880	50.400
1830	43.230	56.130	51.400
				1840	49.090	61.030	52.600
1850	55.430	65.930	54.100
				1860	62.290	71.160	55.000

1870	69.680	76.790	56.500
				1880	77.650	82.550	57.700
1890	86.210	88.630	58.600
				1900	95.390	95.110	59.300
1910	105.240	101.910	59.900
				1920	115.790	109.180	60.500

JP-1 reservoir level ~ storage-capacity curve see Fig. 2, JP-1 reservoir level ~ curve of areas see the long curve in Fig. 3, JP-1 reservoir level ~ storehouse see Fig. 4, according to the form parameter such as water level, storage capacity, storehouse length, area of reservoir, set up model of reservoir as shown in Figure 5.

Wherein: y---water level;

A (y)---the reservoir surface area at y place.

2) water body heat conduction model is set up:

For a water body infinitesimal (see Fig. 5).The thickness that equivalent water level is corresponding is the water body volume of dy is A (y) × dy.

The available following formula of heat motion of this water body represents:

1. Vertical dimension: in the unit interval is Q by the flow entered below _y, the heat brought into is c ρ Q _yt; Be Q by the flow flowed out above _y+ dQ _y, the heat taken away is therefore raw tape enters heat and is in the unit interval

2. within: unit interval, warehouse-in brings heat into level is c ρ _iq _it _idy; It is c ρ q that heat is taken away in outbound ₀tdy; Therefore clean surplus heat is in the unit interval: Q ₂=(c ρ _iq _it _i-c ρ q ₀t) dy.

3. by shortwave radiation heat: from radiation heat R (y) A (y) of leaving away below; From the radiation heat entered above $R\left(y\right)A\left(y\right)+\frac{\∂}{\∂y}\left(RA\right)dy;$ The net radiation heat stayed is: $Q_3=\frac{\∂}{\∂y}\left(RA\right)dy.$

4. by diffusion:

Enter below $-c\mathrm{\ρ}A(D_m+E)\frac{\∂T}{\∂y},$

Top is flowed out $-c\mathrm{\ρ}A(D_m+E)\frac{\∂T}{\∂y}+\frac{\∂}{\∂y}\[-c\mathrm{\ρ}A(D_m+E)\frac{\∂T}{\∂y}\]dy$

Net inflow is $Q_4=\frac{\∂}{\∂y}\[c\mathrm{\ρ}A(D_m+E)\frac{\∂T}{\∂y}\]dy.$

5. water body heats up and absorbs heat: $Q_5=-c\mathrm{\ρ}\frac{\∂T}{\∂\mathrm{\τ}}Ady.$

In various above:

C---specific heat of water;

The density of ρ---water;

The temperature of T---water;

Q _i---the flow of warehouse-in current unit height;

T _i---the temperature of warehouse-in current;

Q ₀---the flow of outbound current unit height;

ρ _i---the density of warehouse-in current;

The shortwave radiation heat at R (y)---height y place;

K---photothermal attenuation coefficient;

A (y)---the reservoir surface area at y place;

D _m---water diffusion coefficient;

The turbulence diffusion coefficient of E---water.

From heat balance:

Q ₁+Q ₂+Q ₃+Q ₄+Q ₅＝0(1-1)

Respectively by Q ₁, Q ₂, Q ₃, Q ₄and Q ₅expression formula bring formula (1-1) into and then have:

$\begin{array}{c}\frac{\∂T}{\∂\mathrm{\τ}}+\frac{1}{c\mathrm{\ρ}A}\frac{\∂}{\∂y}\left({\mathrm{c\ρQ}}_{y}T\right)=\frac{1}{c\mathrm{\ρ}A}\frac{\∂}{\∂y}\[c\mathrm{\ρ}A(D_m+E)\frac{\∂T}{\∂y}\]+\\ +\frac{{\mathrm{\ρ}}_i{q}_i{T}_i}{\mathrm{\ρ}A}-\frac{q_{0}T}{A}+\frac{1}{c\mathrm{\ρ}A}\frac{\∂}{\∂y}\left(RA\right)\end{array}---(1-2)$

3) heat transfer model of warehouse-in, outbound water body is set up:

The density of water body is relevant to water temperature.The density of aqueous water is maximum 4 DEG C time.When water temperature is higher than 4 DEG C, the density of water increases along with the reduction of temperature; In the scope of 0 ~ 4 DEG C, the density of water reduces along with the reduction of temperature, until freezing point.This characteristic of water just, the thermograde that 4 DEG C sink under water, water body are formed from reservoir surface to bottom of the reservior portion from high to low suppresses convection current, even severe cold area, the bottom of reservoir is at least 4 DEG C of stable water layers, can spend severe cold in comfort for hydrobiont.

In the operation of reservoir, warehouse-in current enter the water layer of reservoir equal densities (water temperature) by density (water temperature); And when power station or intake structure diversion, then by corresponding intake height, the water body of the corresponding water temperature layer of carrying-off.The height of water temperature layer, the reservoir shape by the volume and corresponding elevation that go out to put in storage water body is tried to achieve.Therefore, for identical enter the water body of (going out) storehouse volume, in the water body heat interchange of difformity reservoir, the distribution of water temperature change of formation is different.

For formula (1-2), specific heat of water C and water diffusion coefficient D _mfor constant

Q _y＝Q _y(y,τ)，T＝T(y,τ)

And have ρ=ρ ₀+ aT ²+ bT (1-3)

Wherein: the density p of water when 0 DEG C ₀=999.87, a=-0.0067, b=0.07

$\frac{\∂Q_y}{\∂y}=q_i-q_0---(1-4)$

ν(y)＝Q _y/A(1-5)

At y _sthe heat radiation at-y place is:

$R=(1-\mathrm{\β}){\mathrm{\φ}}_0{e}^{-k(y_s-y)}---(1-6)$

In formula: φ ₀for the sun power that the water surface absorbs; β be less than 1.0 coefficient.The turbulence diffusion coefficient E=E (y, τ) of water.

After (1-2) formula of arrangement:

$\begin{array}{c}\frac{\∂T}{\∂\mathrm{\τ}}=\{A\[(D_m+E)\frac{\∂A}{\∂y}+A\frac{\∂E}{\∂y}\]-\frac{v(3{\mathrm{aT}}^2+2bT+{\mathrm{\ρ}}_0)}{{\mathrm{\ρ}}_0+{\mathrm{aT}}^2+bT}\}\frac{\∂T}{\∂y}\\ +\frac{(D_m+E)(2aT+b)}{{\mathrm{\ρ}}_0+{\mathrm{aT}}^2+bT}{\left(\frac{\∂T}{\∂y}\right)}^2+(D_m+E)\frac{{\∂}^{2}T}{\∂y^2}+\frac{q_i}{A}(T_i-T)+\\ +\frac{(1-\mathrm{\β}){\mathrm{\φ}}_0}{c({\mathrm{\ρ}}_0+{\mathrm{aT}}^2+bT)}{e}^{-k(y_s-y)}(\frac{1}{A}\frac{\∂A}{\∂y}+k)\end{array}---(1-7)$

(1-7) is simplified, can following formula be obtained:

$\frac{\∂T}{\∂\mathrm{\τ}}=\{\frac{1}{A}\[(D_m+E)\frac{\∂A}{\∂y}+A\frac{\∂E}{\∂y}\]-v\}\frac{\∂T}{\∂y}+(D_m+E)\frac{{\∂}^{2}T}{\∂y^2}+$

$+\frac{q_i}{A}(T_i-T)+\frac{0.001(1-\mathrm{\β}){\mathrm{\φ}}_0}{c}{e}^{-k(y_s-y)}(\frac{1}{A}\frac{\∂A}{\∂y}+k)---(1-8)$

The basic model that formula (1-8) is hydroelectric station operation phase reservoir water body exchange heat, that is the governing equation of water temperature of reservoir numerical analysis.After the starting condition of given specific reservoir, boundary condition and reservoir service condition, can solve.

For JP-1 reservoir.Its reservoir area meteorological model condition, reservoir service condition are as follows:

Table 2JP-1 hinge flood discharge and diversion facility property list

Table 3JP-1 reservoir goes out reservoir inflow and ruuning situation table

Table 4JP-1 reservoir assignment of traffic table (m ³/ s)

Period	Generator is quoted	Table hole is quoted	Deep hole is quoted	Flood discharging tunnel is quoted	Storage outflow
						January	661				661

February	698				698
						March	753				753
April	832				832
						May	916				916
June	1455				1455
						July	2076	754			1854
August	2076	524			2169
						September	2076	444			2177
October	1591				1591
						November	810				810
Dec	640				640

Note: overflow weir crest level power station irrigating gate center line elevation

Table hydrometric station, 5JP-1 reservoir area meteorological element statistical form

In the present embodiment, set up the mathematical model of water temperature of reservoir analysis according to heat-conduction principle after, in solution approach, implicit difference scheme can be adopted.For the Solve problems of one dimension reservoir water temperature model, be summed up as and solve a tri-diagonal matrix equation, after carrying out LU decomposition, try to achieve asked a question solution.

Governing equation (1-8) using differential method solution it.In order to avoid the restriction of explicit scheme difference equation stable condition, this method adopts implicit difference solution.

Order $\left.\begin{array}{c}\frac{\∂T}{\∂\mathrm{\τ}}=\frac{T(y,\mathrm{\τ})-T(y,\mathrm{\τ}-\mathrm{\Δ}\mathrm{\τ})}{\mathrm{\Δ}\mathrm{\τ}}\\ \frac{\∂T}{\∂y}=\frac{T(y+\mathrm{\Δ}y,\mathrm{\τ})-T(y-\mathrm{\Δ}y,\mathrm{\τ})}{\text{2}\mathrm{\Δ}y}\\ \frac{{\∂}^{2}T}{\∂y^2}=\frac{T(y+\mathrm{\Δ}y,\mathrm{\τ})+T(y-\mathrm{\Δ}y,\mathrm{\τ})-2T(y,\mathrm{\τ})}{{\mathrm{\Δy}}^2}\\ \frac{\∂A}{\∂y}=\frac{A(y+\mathrm{\Δ}y)-A\left(y\right)}{\mathrm{\Δ}y}\\ \frac{\∂E}{\∂y}=\frac{E(y+\mathrm{\Δ}y,\mathrm{\τ})-E(y,\mathrm{\τ})}{\mathrm{\Δ}y}\end{array}\right\}---(1-9)$

Substitute into (1-8) arrangements must:

$\begin{array}{c}\frac{\mathrm{\Δ}\mathrm{\τ}}{{\mathrm{\Δy}}^2}\{(D_m+E)+\frac{1}{2}\[(D_m+E)\frac{\mathrm{\Δ}A}{A\left(y\right)}+\mathrm{\Δ}E-\mathrm{\Δ}yv(y,\mathrm{\τ})\]\}T(y+\mathrm{\Δ}y,\mathrm{\τ})+\\ +\{-\frac{2(D_m+E)\mathrm{\Δ}\mathrm{\τ}}{{\mathrm{\Δy}}^2}-\frac{{\mathrm{\Δ\τq}}_i(y,\mathrm{\τ})}{A\left(y\right)}-1\}T(y,\mathrm{\τ})+\\ +\{\frac{\mathrm{\Δ}\mathrm{\τ}}{2{\mathrm{\Δy}}^2}\[2(D_m+E)-(D_m+E)\frac{\mathrm{\Δ}A}{A\left(y\right)}-\mathrm{\Δ}E+\mathrm{\Δ}yv(y,\mathrm{\τ})\]\}T(y-\mathrm{\Δ}y,\mathrm{\τ})\\ =-\frac{{\mathrm{\Δ\τq}}_i(y,\mathrm{\τ})}{A\left(y\right)}{T}_i(y,\mathrm{\τ})-\frac{0.001\mathrm{\Δ}\mathrm{\τ}(1-\mathrm{\β}){\mathrm{\φ}}_0}{c\mathrm{\Δ}y}{e}^{-k(y_s-y)}\[\frac{\mathrm{\Δ}A}{A\left(y\right)}+k\mathrm{\Δ}y\]-T(y,\mathrm{\τ}-\mathrm{\Δ}\mathrm{\τ})\end{array}---(1-10)$

In formula:

$\left.\begin{array}{c}\mathrm{\Δ}A=A(y+\mathrm{\Δ}y)-A\left(y\right)\\ \mathrm{\Δ}E=E(y+\mathrm{\Δ}y,\mathrm{\τ})-E(y,\mathrm{\τ})\end{array}\right\}---(1-11)$

It is more than the equation that arbitrary y point is set up.If so water body is divided into n layer, then have 1,2 ..., n, n+1 common (n+1) individual point.If make Δ y=y _s/ n, then y _s=n Δ y, so to general some j point, then have:

$\begin{array}{c}\frac{\mathrm{\Δ}\mathrm{\τ}}{{\mathrm{\Δy}}^2}\{\[D_m+E\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\right)\]+\[\[D_m+E\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\right)\]\frac{A\left(j\mathrm{\Δ}y\right)-A\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y\right)}{A\left(j\mathrm{\Δ}y\right)}+\\ +\[E(j\mathrm{\Δ}y,\mathrm{\τ})-E\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\right)\]-\mathrm{\Δ}yv\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\right)\]\}T(j\mathrm{\Δ}y,\mathrm{\τ})+\\ +\left\{-\frac{2\mathrm{\Δ}\mathrm{\τ}}{{\mathrm{\Δy}}^2}\[D_m+E\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\right)\]-\frac{{\mathrm{\Δ\τq}}_i\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\right)}{A\left(j\mathrm{\Δ}y\right)}-1\right\}T\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\right)+\\ +\{\frac{\mathrm{\Δ}\mathrm{\τ}}{2{\mathrm{\Δy}}^2}\[2\left(D_m+E\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\right)\right)-\left(D_m+E\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\right)\right)\frac{A\left(j\mathrm{\Δ}y\right)-A\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y\right)}{A\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y\right)}-\\ -\left(E\left(j\mathrm{\Δ}y,\mathrm{\τ}\right)-E\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\right)\right)+\mathrm{\Δ}yv\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\right)\]\}T\left(\stackrel{\‾}{j-2}\mathrm{\Δ}y,\mathrm{\τ}\right)\\ =-\frac{{\mathrm{\Δ\τq}}_i\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\right)}{A\left(j\mathrm{\Δ}y\right)}{T}_i\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\right)-\frac{0.001\mathrm{\Δ}\mathrm{\τ}\left(1-\mathrm{\β}\right){\mathrm{\φ}}_0}{c\mathrm{\Δ}y}{e}^{-k\left(n-j+1\right)\mathrm{\Δ}y}.\\ \·\[k\mathrm{\Δ}y+\frac{A\left(j\mathrm{\Δ}y\right)-A\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y\right)}{A\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y\right)}\]-T\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}-\mathrm{\Δ}\mathrm{\τ}\right)\end{array}---(1-12)$

In (1-12) formula, equal sign right-hand member is given value, equal sign left end ... in be known, only have temperature T (j Δ y, τ), with for the unknown.So to every bit j (j=1 ..., n+1), an above-mentioned equation can be set up, total n+1 equation, just can solve n+1 temperature T (j Δ y, τ) (j=1 ..., n+1).Being write as matrix form is:

[K]{T}＝{F}(1-13)

In above formula matrix equation:

$\left\{T\right\}={[T_1(0,\mathrm{\τ}),T_2(\mathrm{\Δy},\mathrm{\τ})...T_j(\stackrel{\‾}{j-1}\mathrm{\Δy},\mathrm{\τ})...T_{n+1}(\mathrm{n\Δy},\mathrm{\τ})]}^T---(1-14)$

{F}＝[f ₁f ₂…f _n+1] ^T(1-15)

$\[K\]=\left[\begin{array}{cccccc}{K}_{11}& K_{12}& & & & 0\\ K_{21}& K_{22}& K_{23}& & & \\ & K_{32}& K_{33}& K_{34}& & \\ & & K_{j,j-1}& K_{jj}& K_{j,j+1}& \\ & & & K_{n,n-1}& K_{nn}& K_{n,n+1}\\ & & & & K_{n+1,n}& K_{n+1,n+1}\end{array}\right]---(1-16)$

Each element in formula (1-15), (1-16), is tried to achieve by (1-12) formula, for:

$f_j=-\frac{{\mathrm{\Δ\τq}}_i(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ})}{A\left(j\mathrm{\Δ}y\right)}{T}_i(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ})-\frac{0.001\mathrm{\Δ}\mathrm{\τ}(1-\mathrm{\β}){\mathrm{\φ}}_0}{c\mathrm{\Δ}y}{e}^{-k(n-j+1)\mathrm{\Δ}y}\·$

$\·\[k\mathrm{\Δ}y+\frac{\text{A(j}\mathrm{\Δ}\text{y)-A}\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y\right)}{\text{A}\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y\right)}\]-T(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}-\mathrm{\Δ}\mathrm{\τ}),(j=2,...,n)\text{---}(1-17)$

$K_{jj}=-\frac{2\mathrm{\Δ}\mathrm{\τ}}{{\mathrm{\Δy}}^2}\[D_m+E(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ})\]+\frac{{\mathrm{\Δ\τq}}_i(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ})}{A\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y\right)}-1---(1-18)$

$\begin{array}{c}{K}_{j,j-1}=\frac{\mathrm{\Δ}\mathrm{\τ}}{2{\mathrm{\Δy}}^2}\[2(D_m+E(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\left)\right)-(D_m+E(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\left)\right)\·\\ \·\frac{A\left(j\mathrm{\Δ}y\right)-A\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y\right)}{A\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y\right)}-\left(E\right(j\mathrm{\Δ}y,\mathrm{\τ}\left)\right)-E\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\right))+\mathrm{\Δ}yv(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ})\]\end{array}---(1-19)$

$\begin{array}{c}{K}_{j,j+1}=\frac{\mathrm{\Δ}\mathrm{\τ}}{{\mathrm{\Δy}}^2}\{D_m+E(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ})+\frac{1}{2}\[(D_m+E(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\left)\right)\·\\ \·\frac{A\left(j\mathrm{\Δ}y\right)-A\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y\right)}{A\left(j\mathrm{\Δ}y\right)}+\left(E\right(j\mathrm{\Δ}y,\mathrm{\τ}\left)\right)-E\left(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ}\right))-\\ -\mathrm{\Δ}yv(\stackrel{\‾}{j-1}\mathrm{\Δ}y,\mathrm{\τ})\]\}\end{array}---(1-20)$

In above three formulas, j=2 ..., n.

About f ₁, f _n+1and K ₁₁, K ₁₂, K _{n+1, n}, K _{n+1, n+1}etc. belonging to boundary condition treatment, be omitted herein.

For the border of bottom of the reservior, narrowly it is alternating temperature border, but considers bottom of the reservior certain depth water layer, and its temperature is generally more stable, and calculate for simplifying, be similar to and can be taken as adiabatic condition process, error is very little.So can to downward-extension one lattice, this temperature should equal y=Δ y place temperature.Therefore have:

$\begin{array}{c}{f}_1=\frac{-{\mathrm{\Δ\τq}}_i(0,\mathrm{\τ})}{A\left(0\right)}{T}_i(0,\mathrm{\τ})-\frac{0.001\mathrm{\Δ}\mathrm{\τ}(1-\mathrm{\β}){\mathrm{\φ}}_0}{c\mathrm{\Δ}y}{e}^{-{\mathrm{ky}}_s}\·\\ \·\[\frac{A\left(\mathrm{\Δ}y\right)-A\left(0\right)}{A\left(0\right)}+k\mathrm{\Δ}y\]-T(0,\mathrm{\τ}-\mathrm{\Δ}\mathrm{\τ})\end{array}---(1-21)$

$K_{11}=-\frac{2\mathrm{\Δ}\mathrm{\τ}}{{\mathrm{\Δy}}^2}\[D_m+E(0,\mathrm{\τ})\]-\frac{{\mathrm{\Δ\τq}}_i(0,\mathrm{\τ})}{A\left(0\right)}-1---(1-22)$

$K_{12}=\frac{2\mathrm{\Δ}\mathrm{\τ}}{{\mathrm{\Δy}}^2}\[D_m+E(0,\mathrm{\τ})\]---(1-23)$

For the boundary condition of the water surface, can be write as:

$\mathrm{\ρ}c\[D_m+E(y,\mathrm{\τ})\]\frac{\∂T_s}{\∂y}=-{\mathrm{\β\φ}}_0-{\mathrm{\φ}}_a+{\mathrm{\φ}}_b+{\mathrm{\φ}}_e+{\mathrm{\φ}}_c---(1-24)$

In formula: β φ ₀---the sun power that water surface absorbs; φ _a---the atmosphere radiation energy that the water surface absorbs; φ _b---the water surface returns the radiation energy of air; φ _e---the latent heat that evaporation from water surface absorbs; φ _c---water surface convection losses heat.Its formula is as follows:

$\mathrm{\β}{\mathrm{\φ}}_0=0.94\mathrm{\β}{\mathrm{\φ}}_{\mathrm{sc}}(1.0-0.65c_c^2)$ (1-25)

Wherein: φ _sc---solar radiant energy during fine day, substantially relevant with the latitude of position, if can search by [4] without during field data.Cc is the cloud-capped percentage of sky.

${\mathrm{\φ}}_a=0.97\mathrm{\ϵ\σ}{(T_a+273)}^4(1.0+0.17c_c^2)$ (1-26)

In formula: ε---air average radiation energy;

$\mathrm{\ϵ}=1.0-0.26/\exp \{7.77\×{10}^{-5}\left(T_a^2\right)\}$ (1-27)

Here: σ---be Stefen-Boltzmann constant,

σ＝4.9×10 ^-3J/m ²·d·°K ⁴＝0.001170×10 ^-3Kcal/m· ²d·°K ⁴；

T _a---with DEG C temperature represented.

φ _b＝0.97σ(T _s+273) ⁴＝1.135×10 ^-6(T _s+273) ⁴Kcal/m ²·d(1-28)

φ _e＝(2493-2.26T _s)×10 ³E(Kcal/m ²·d)(1-29)

In formula: E=Δ h × 10 ⁴/ 4186.8 × 12 × 30.4, Δ h in mm, is evaporation capacity.

φ _c＝269.1ρ(T _s-T _z)(0.000308+0.000185W _z)(1-30)

Tz, Wz are respectively temperature and the wind speed at high z place on the water surface, the value at high 15cm (6 inch) place on general water intaking face.Above all formulas are substituted into (1-29) and arrange it, must water surface temperature T be determined _s=T (y _s, τ) formula:

T _s ⁴+1092T _s ³+447174T _s ²+(81385668-1991364878E+2.371134×10 ¹¹α ₁)T _s+5554571841+2.1966693×10 ¹²E-2.371134×10 ¹¹α ₁T _z-α ₂-α ₃＝0(1-31)

In formula:

$\left.\begin{array}{c}{\mathrm{\α}}_1=0.000308+0.000185W_z\\ {\mathrm{\α}}_2=828266.8077{\mathrm{\β\φ}}_{sc}(1-0.65c_c^2)\\ {\mathrm{\α}}_3=\mathrm{\ϵ}{(T_a+273)}^4(1+0.17c_c^2)\end{array}\right\}---(1-32)$

When deducing (1-31), assuming that T (y _s, τ) and=T (y _s-Δ y, τ).When Δ y value is little, consider that the water surface has blending effect, this hypothesis allows (in this method calculating, Δ y≤0.5 meter).

Formula (1-31) is standard unary biquadratic equation, and general available equation, dichotomy or Secant Method solve.This method adopts dichotomy to solve.So T (y _s, τ) and for known, thus can obtain:

$K_{nn}=-\frac{2\mathrm{\Δ}\mathrm{\τ}}{{\mathrm{\Δy}}^2}\[D_m+E(y_s-\mathrm{\Δ}y,\mathrm{\τ})\]+\frac{{\mathrm{\Δ\τq}}_i(y_s-\mathrm{\Δ}y,\mathrm{\τ})}{A(y_s-\mathrm{\Δ}y)}-1\text{---}(1-33)$

$\begin{array}{c}{K}_{n,n-1}=\frac{\mathrm{\Δ\τ}}{2\mathrm{\Δ}{y}^2}[2(D_m+E(y_s-\mathrm{\Δy},\mathrm{\τ})-(D_m+E(y_s-\mathrm{\Δy},\mathrm{\τ}))\\ \frac{A\left(y_s\right)-A(y_s-\mathrm{\Δy})}{A(y_s-\mathrm{\Δy})}-(E(y_s,\mathrm{\τ})-E(y_s-\mathrm{\Δy},\mathrm{\τ}))+\\ +\mathrm{\Δyv}(y_s-\mathrm{\Δy},\mathrm{\τ})\end{array},---(1-34)$

$\begin{array}{c}{f}_n=-\frac{{\mathrm{\Δ\τq}}_i(y_s-\mathrm{\Δ}y,\mathrm{\τ})}{A(y_s-\mathrm{\Δ}y)}{T}_i(y_s-\mathrm{\Δ}y,\mathrm{\τ})-\\ -\frac{0.001\mathrm{\Δ}\mathrm{\τ}(1-\mathrm{\β}){\mathrm{\φ}}_0}{c\mathrm{\Δ}y}{e}^{-k\mathrm{\Δ}y}\[k\mathrm{\Δ}y+\frac{A\left(y_s\right)-A(y_s-\mathrm{\Δ}y)}{A(y_s-\mathrm{\Δ}y)}\]-\\ -T(y_s-\mathrm{\Δ}y,\mathrm{\τ}-\mathrm{\Δ}\mathrm{\τ})-\frac{\mathrm{\Δ}\mathrm{\τ}}{{\mathrm{\Δy}}^2}\{\[D_m+E(y_s-\mathrm{\Δ}y,\mathrm{\τ})\]+\\ +\frac{1}{2}\[(D_m+E(y_s-\mathrm{\Δ}y,\mathrm{\τ}\left)\right)\frac{A\left(y_s\right)-A(y_s-\mathrm{\Δ}y)}{A(y_s-\mathrm{\Δ}y)}+\\ +\left(E\right(y_s,\mathrm{\τ})-E(y_s-\mathrm{\Δ}y,\mathrm{\τ}\left)\right)-\mathrm{\Δ}yv\left(y_s-\mathrm{\Δ}y,\mathrm{\τ}\right)\]\}T(y_s,\mathrm{\τ})\end{array}---(1-35)$

The matrix of (1-16) formula reduces a line one row like this:

$\[K\]=\left[\begin{array}{cccccc}{K}_{11}& K_{12}& & & & 0\\ K_{21}& K_{22}& K_{23}& & & \\ & K_{32}& K_{33}& K_{34}& & \\ & & K_{j,j-1}& K_{jj}& K_{j,j+1}& \\ & & & K_{n-1,n-2}& K_{n-1,n-1}& K_{n-1,n}\\ 0& & & & K_{n,n-1}& K_{nn}\end{array}\right]---(1-36)$

Above formula is three diagonal angle band matrix; For saving memory capacity during calculating, only can store the nonzero element on three diagonal line, being about to [K] matrix [V] matrix and replacing as follows:

$\[V\]=\left[\begin{array}{ccc}{Q}_1& V_{12}& V_{13}\\ V_{21}& V_{22}& V_{23}\\ V_{31}& V_{32}& V_{33}\\ \·& \·& \·\\ \·& \·& \·\\ \·& \·& \·\\ V_{i\text{1}}& V_{i2}& V_{i3}\\ \·& \·& \·\\ \·& \·& \·\\ \·& \·& \·\\ V_{n,1}& V_{n,2}& Q_2\end{array}\right]---(1-37)$

Q1, Q2 in formula are Arbitrary Digit.Element Kij and element V in [V] in [K] _{i ' j '}pass be:

Present m=3, so j '=j-i+2.Such as K11 is V12 in [V]; Knn is Vn in [V], 2.

About the banded solution of equations method in three diagonal angles, first can make [L] [U] to [K] and decompose:

[K]＝[L][U](1-39)

Wherein:

$\[L\]=\left[\begin{array}{cccc}{\mathrm{\β}}_1& & & 0\\ K_{21}& {\mathrm{\β}}_2& & \\ & K_{32}& {\mathrm{\β}}_3& \\ & & K_{n,n-1}& {\mathrm{\β}}_n\end{array}\right]\text{---}(1-40)$

$\[U\]=\left[\begin{array}{ccccc}1& {\mathrm{\δ}}_1& & & 0\\ & 1& {\mathrm{\δ}}_2& & \\ & & 1& {\mathrm{\δ}}_3& \\ & & & 1& {\mathrm{\δ}}_{n-1}\\ 0& & & & 1\end{array}\right]\text{---}(1-41)$

Here:

$\left.\begin{array}{c}{\mathrm{\β}}_1=K_{11}\\ {\mathrm{\δ}}_k=K_{k,k+1}/{\mathrm{\β}}_k,k=1,2,...,n-1\\ {\mathrm{\β}}_{k+1}=K_{k+1,k+1}-K_{k,k+1}{\mathrm{\δ}}_k\end{array}\right\}---(1-42)$

So (1-13) equivalent equation of formula is:

$\left.\begin{array}{c}\[L\]\left\{Z\right\}=\left\{F\right\}\\ \[U\]\left\{T\right\}=\left\{Z\right\}\end{array}\right\}---(1-43)$

Its formula is as follows:

$\left.\begin{array}{cc}{z}_1=f_1/K_{11}& \\ z_k=(f_k-K_{k,k-1}z{k-1}_{})/{\mathrm{\β}}_k& (k=2,...,n)\end{array}\right\}---(1-44)$

$\left.\begin{array}{cc}{T}_n=z_n& \\ T_k=z_k-{\mathrm{\δ}}_k{T}_{k+1},& (k=n-1,...,1)\end{array}\right\}---(1-45)$

Store because [K] matrix takes [V] matrix, corresponding computing formula is as follows:

$\left.\begin{array}{c}{\mathrm{\β}}_1=K_{11}=V_{12}\\ {\mathrm{\δ}}_k=K_{k,k+1}/{\mathrm{\β}}_k=V_{k,3}/{\mathrm{\β}}_k\\ {\mathrm{\β}}_{k+1}=K_{k+1,k+1}-K_{k,k-1}{\mathrm{\δ}}_k=V_{k+1,2}-V_{k,1}{\mathrm{\δ}}_k\\ (k=1,2,...,n-1)\end{array}\right\}---(1-46)$

$\left.\begin{array}{c}{z}_1=f_1/K_{11}=f_1/V_{12}\\ z_k=\left(f_k-K_{k,k-1}{z}_{k-1}\right)/{\mathrm{\β}}_k=\left(f_k-V_{k,1}{z}_{k-1}\right)/{\mathrm{\β}}_k\}\\ \left(k=2,...,n\right)\end{array}\right\}---(1-47)$

$\left.\begin{array}{cc}{T}_n=z_n& \\ T_k=z_k-{\mathrm{\δ}}_k{T}_{k+1},& \left(k=n-1,...,1\right)\end{array}\right\}---(1-48)$

So, { after F}, calculate { β }, { δ } by formula (1-46) at known [V], calculated by formula (1-47) that { z} by formula (1-48), can calculate the water temperature { T} of Different periods reservoir different depth in year.

In the present embodiment, calculate the process of centering to Railway Project to comprise:

(1) simulation of water volume flow rate ν:

$v(y,\mathrm{\τ})=\frac{Q_y(y,\mathrm{\τ})}{A\left(y\right)}=\frac{1}{A\left(y\right)}\underset{0}{\overset{y}{\∫}}\[q_i(\mathrm{\η},\mathrm{\τ})-q_0(\mathrm{\η},\mathrm{\τ})\]d\mathrm{\η}---(1-49)$

(2) water level y _scalculating:

If a certain initial time τ ₀, reservoir level is y ₀, storehouse water capacity is V ₀, to moment τ ₁=τ ₀water level during+Δ τ is y _s, storehouse water capacity is within this Δ τ time interval, reservoir inflow is Q _i(τ ₁), storage outflow is Q ₀(τ ₁), then have:

V _τ1＝V ₀+[Q _i(τ _i)-Q ₀(τ ₁)]Δτ(1-50)

Can by water level storage-capacity curve, according to value checks in y _s.

(3) initial temperature: generally can fix on spring, water temperature of reservoir is evenly T ₀.

(4) current are put in storage:

The velocity flow profile that can be similar to hypothesis warehouse-in current is normal distribution:

$u_i(y,\mathrm{\τ})=u_{im}\left(\mathrm{\τ}\right)e^{-{(y-y_i)}^2/2{\mathrm{\σ}}_i^2}---(1-51)$

If warehouse-in total flow is Q _i(τ), then have

$Q_i\left(\mathrm{\τ}\right)=u_m\left(\mathrm{\τ}\right){\∫}_0^{y_s}B\left(y\right)e^{-{(y-y_i)}^2/2{\mathrm{\σ}}_i^2}dy---(1-52)$

In above formula: σ _ifor standard deviation; u _im(τ) be the Peak Flow Rate of reservoir current by normal distribution during τ, the position of yi Peak Flow Rate for this reason.

B(y)＝A(y)/L(y)(1-53)

In formula: the reservoir average length that L (y) is y place.

(1-53) is substituted into (1-52), and is write as the form of Cumulative sum, can obtain:

$u_{im}\left(\mathrm{\τ}\right)=Q_i\left(\mathrm{\τ}\right)/\underset{j=1}{\overset{n}{\Σ}}\frac{\mathrm{\Δ}yA\left(j\mathrm{\Δ}y\right)}{L\left(j\mathrm{\Δ}y\right)}{e}^{-{(j\mathrm{\Δ}y-y_i)}^2/2{\mathrm{\σ}}_{\text{i}}^2}---(1-54)$

Suppose that in natural river course, even water temperature is Ti (due to river action of turbulent flow), after warehouse-in, the center elevation of current is the elevation (putting aside that warehouse-in current are entering the blending effect at storehouse place) of storehouse coolant-temperature gage Ti.If mixed water current, then before first checking the mixed water being less than certain particle that dam can be arrived, and enter the equal densities layer of reservoir by its density.The thickness of warehouse-in current in storehouse is δ _i, be similar to and be taken as:

${\mathrm{\δ}}_i=4.8{\left(\frac{q_{ii}^2}{g\mathrm{\ϵ}}\right)}^{\text{1/4}}---(1-55)$

q _ii＝Q _i(τ)/B _i(1-56)

B _i＝A(y _i)/L(y _i)(1-57)

Observe according in flume test and general reservoir, regulation:

${\mathrm{\σ}}_i=\frac{0.5{\mathrm{\δ}}_i}{1.96}---(1-58)$

In (1-55): g is acceleration of gravity;

$\mathrm{\ϵ}=\frac{1}{\mathrm{\ρ}}\·\frac{\∂\mathrm{\ρ}}{\∂y}---(1-59)$

Thus have:

$\mathrm{\ϵ}=\frac{1}{\mathrm{\ρ}}(2aT+b)\frac{\∂T}{\∂y}---(1-60)$

Write as difference form to have:

${\mathrm{\δ}}_i=4.8{\left\{\frac{q_{ii}^2\[{\mathrm{aT}}^2(y_i,\mathrm{\τ})+bT(y_i,\mathrm{\τ})+{\mathrm{\ρ}}_0\]\·2\mathrm{\Δ}y}{g\[2aT(y_i,\mathrm{\τ})+b\]\·\[T(y_i+\mathrm{\Δ}y,\mathrm{\τ})-T(y_i-\mathrm{\Δ}y,\mathrm{\τ})\]}\right\}}^{1/4}---(1-61)$

U can be tried to achieve like this by above-mentioned formula _i(y, τ), thus q _i(y, τ) is known.For outbound current, can be imitated this and try to achieve.

The present embodiment, for JP-1 reservoir, utilizes above-mentioned Forecasting Methodology, and gained numerical results is see table 6, table 7, Fig. 6, Fig. 7, Fig. 8, Fig. 9.

Table 6JP-1 water temperature of reservoir distribution (not considering front bottom of the reservior, dam slag, unit DEG C)

Table 7JP-1 water temperature of reservoir distribution (considering front bottom of the reservior, dam slag, unit DEG C)

Month	January	February	March	April	May	June	July	August	September	October	November	Dec	Every year
														Temperature	10.5	13.6	17.4	20.5	21.7	21.9	21.2	21.4	19.2	17.1	12.8	9.1	17.2
Water level (m)	1868	1854	1836	1815	1802	1800	1846	1880	1880	1880	1880	1877	1851.5
														1880								22.7	20.7	18.6	15.1	11.5
1875								22.6	20.2	18.6	15.0	11.5
														1870	12.5							22.6	19.9	18.5	14.9	11.5
1865	12.0							22.6	19.3	18.4	14.8	11.5
														1860	11.8							22.3	18.9	18.2	14.7	11.5
1855	11.7	14.9						22.2	18.6	17.9	14.7	11.5
														1850	11.6	14.7						21.8	18.3	17.8	14.6	11.5	18.6
1845	11.6	14.2					22.8	21.5	17.8	17.5	14.5	11.5	18.2
														1840	11.5	13.9					22.6	21.2	17.7	17.3	14.4	11.5	18.0
1835	11.5	13.4	19.0				22.4	21.0	17.5	17.0	14.4	11.5	17.9
														1830	11.5	13.0	18.8				22.0	20.6	17.4	16.9	14.3	11.5	17.8
1825	11.5	12.5	18.2				21.5	20.0	17.3	16.8	14.3	11.5	17.6
														1820	11.4	12.2	17.8				20.8	19.8	17.2	16.5	14.2	11.5	17.4
1815	11.5	12.0	16.8	21.0			20.2	19.5	17.1	16.5	14.1	11.5	17.2
														1810	11.5	11.8	15.6	20.9			19.8	19.3	17.1	16.4	14.1	11.5	17.0

1805	11.5	11.7	14.8	20.3			19.5	19.3	17.1	16.4	14.0	11.5	16.8
														1800	11.5	11.7	14.0	19.6	22.3	23.0	19.4	19.3	17.1	16.2	13.9	11.5	16.6
1795	11.5	11.7	13.0	18.8	22.0	22.4	19.4	19.3	17.0	16.1	13.9	11.5	16.4
														1790	11.5	11.7	12.5	17.6	21.0	21.5	19.3	19.3	16.9	16.1	13.6	11.5	16.0
1785	11.5	11.7	12.1	15.6	18.8	20.7	19.2	19.3	16.8	16.1	13.2	11.5	15.5
														1783	11.5	11.6	12.0	15.1	18.2	20.4	19.1	19.3	16.8	16.1	13.1	11.4	15.4
1780	11.4	11.6	11.7	14.2	16.8	19.8	19.0	19.1	16.8	16.1	12.9	11.4	15.1
														1775	11.3	11.4	11.5	13.73	16.0	18.62	18.8	18.7	16.8	16.1	12.6	11.4	14.7
1770	11.0	11.4	11.4	13.5	15.3	17.87	18.3	18.6	16.8	16.1	12.3	11.3	14.5
														1765	10.8	11.1	11.3	12.8	15.0	16.8	17.6	18.2	16.8	16.1	12.0	11.1	14.1
1760	10.5	10.7	11.0	12.5	14.2	15.44	16.7	17.7	16.8	16.1	11.6	10.9	13.7
														1755	10.0	10.32	10.5	12.0	13.8	14.51	15.9	17.0	16.6	16.1	11.0	10.5	13.2
1750	9.7	10.2	10.4	11.6	13.2	13.85	15.5	16.5	16.5	16.0	10.8	10.0	12.9
														1745	9.6	10.0	10.2	11.0	12.3	13.1	14.4	16.2	16.1	15.4	10.2	9.9	12.4
1740	9.4	9.4	9.5	10.5	11.8	12.35	13.8	15.7	15.3	15.1	10.1	9.83	11.9
														1735	9.2	9.3	9.4	10	10.9	11.42	13.2	14.9	14.7	14.5	10.0	9.7	11.4
1730	9.0	9.1	9.1	9.6	10.4	11.0	13.0	14.1	14.4	13.9	10.0	9.6	11.1
														1725	9.0	9.2	9.2	9.5	10.0	10.8	12.5	13.5	13.8	13.2	10.0	9.5	10.8
1720	9.1	9.2	9.3	9.4	9.8	10.5	12.0	13.0	13.3	12.8	10.0	9.4	10.6
														1700	9.0	9.3	9.3	9.4	9.5	9.8	11.0	11.6	11.4	11.4	10.0	9.4	10.1
1680	9.0	9.4	9.5	9.5	9.55	9.75	10.6	11.0	10.6	10.6	10.1	9.5	9.9
														1660	9.3	9.6	9.7	9.7	9.8	10.0	10.4	10.6	10.5	10.6	10.5	9.7	10.0
1640	10.2	10.5	10.6	10.6	10.6	11.0	11.0	11.0	11.2	11.2	10.6	10.6	10.8
														1630	10.9	11.1	11.5	11.5	11.5	11.6	11.6	11.6	11.6	11.7	11.5	11.5	11.5
1620	12.2	12.3	12.8	12.8	12.87	13.0	13.0	13.0	13.0	13.0	12.9	12.8	12.8
														1600	14.1	14.5	14.6	14.6	14.6	14.7	14.7	14.7	14.7	14.7	14.6	14.6	14.6
1580	14.0	14.2	14.3	14.3	14.3	14.4	14.4	14.4	14.4	14.4	14.3	14.3	14.3

The Forecasting Methodology that the present embodiment provides can take into full account the influence factor such as shape, the initial retaining situation of reservoir area hydrometeorology condition, Reservoir Operation Conditions and reservoir of reservoir.Thus can for different reservoir feature separately, the distribution of water temperature rule that forecast analysis Reservoir is later, predicts the outcome rationally accurately.

A series of detailed description listed is above only illustrating for feasibility embodiment of the present invention; they are also not used to limit the scope of the invention, all do not depart from the skill of the present invention equivalent implementations done of spirit or change all should be included within protection scope of the present invention.

To those skilled in the art, obviously the invention is not restricted to the details of above-mentioned one exemplary embodiment, and when not deviating from spirit of the present invention or essential characteristic, the present invention can be realized in other specific forms.Therefore, no matter from which point, all should embodiment be regarded as exemplary, and be nonrestrictive, scope of the present invention is limited by claims instead of above-mentioned explanation, and all changes be therefore intended in the implication of the equivalency by dropping on claim and scope are included in the present invention.

Claims (3)

1. a large hydropower station storehouse coolant-temperature gage Numerical Predicting Method, is characterized in that, comprise the steps:

Step 1, gathers the underlying parameter of Dam Site; Described underlying parameter comprises the initial retaining parameter of reservoir form parameter, hydrometeorological parameter, reservoir operational factor and reservoir;

Step 2, the underlying parameter according to gathering sets up numerical analysis model;

Step 3, is defined as low flow year parameter, normal flow year parameter and high flow year parameter according to low flow year, normal flow year and high flow year by the underlying parameter of collection;

Step 4, according to low flow year parameter, normal flow year parameter and high flow year parameter, the water temperature Changing Pattern respectively in numerical simulation year under reservoir laminar flow exchange heat;

Step 5, again respectively for not considering upstream Ba Qian bottom of the reservior slag and alluvial, the alluvial of Ba Qian bottom of the reservior, consideration upstream, considering upstream Ba Qian bottom of the reservior slag and alluvial three kinds of situations, carry out numerical simulation, obtain low value result, median result, the high level result of reservoir runtime upstream Ba Qian storehouse coolant-temperature gage;

Step 6, according to the low value result of reservoir runtime upstream Ba Qian storehouse coolant-temperature gage obtained, median result and high level result, determines the regularity of distribution and the value of runtime water temperature of reservoir, obtains water temperature of reservoir envelop-type and predicts the outcome.

2. a kind of large hydropower station storehouse according to claim 1 coolant-temperature gage Numerical Predicting Method, is characterized in that, described reservoir form parameter comprises relation between reservoir capacity, the reservoir degree of depth, reservoir level ~ storage capacity ~ storehouse length ~ area; Described hydrometeorological parameter comprises temperature, solar radiation, wind speed, cloud amount, evaporation capacity, reservoir inflow, warehouse-in water temperature, river load content, warehouse-in suspended load; Described reservoir operational factor comprises dam safety evaluation mode, irrigating gate position, power station and diversion ability, flash buildings position and discharging capacity, the traffic control situation of reservoir, storage outflow, reservoir level change; Slag situation, upstream cofferdam disposition before ground temperature, first filling temperature, reservoir filling speed, dam when the initial retaining parameter of described reservoir comprises first filling season, first filling.

3. a kind of large hydropower station storehouse according to claim 1 coolant-temperature gage Numerical Predicting Method, it is characterized in that, described step 2 specifically comprises:

Model of reservoir is set up according to reservoir form parameter;

Set up water body heat conduction model;

Set up the heat transfer model of warehouse-in, outbound water body.