Software for Pollutants Transport in Rivers and for Identification of Excessive Pollution Sources
Transport and diffusion equation
The diffusion of polluting substances in rivers is mostly described by the three-dimensional equation of turbulent diffusion of non-conservative substances [1-4]:
is the time-averaged concentration of non-conservative pollutant; t is the time;
is the radius vector; it’s components x,y,z are the spatial coordinates (the axis x is horizontal and its direction coincides with the direction of averaged current of the flow, the axis
is perpendicular to the free surface and is directed downwards; the axis
is directed across the flow);
is the tensor of turbulent diffusion,
is the vector of time-averaged speed of the river flow;
is the coefficient of non-conservatism of pollutants,
is the power of polluting sources.
For the engineering practice very often the two-dimensional or one-dimensional models are considered instead of the three-dimensional one [2,5]. If the non-uniformity of distribution of concentrations of pollutants on the depth of the watercourse is not taken into account, then it is possible to obtain the two-dimensional turbulent diffusion equation. The equation of one-dimensional turbulent diffusion is applied when the distribution of the concentration of the pollutant across the stream is homogeneous. It is also used when average indicators are used to represent pollution across the river .
Initial and boundary conditions
The unknown function
is defined at
is some region with the boundary
(in one-dimensional model).
The additional conditions are specified in the form of
). The boundary conditions in the lower extremity of the river section may be classical or non-classical. The classical conditions (condition of complete mixing) look like
; (condition of full mixing)
non-classical conditions –
. (not local boundary condition)
is the coefficient of self-cleaning of the river on the considered section,
is the length of the section [5,6].
At using of two- or three-dimensional models, the following boundary conditions on the other part of the line or on the surface
(Neumann condition), are set also:
is unit vector of external normal to the border
. In particular, at
, it should be
Algorithms of solving of the initial-boundary-value problems
For solving the considered initial-boundary-value problems one of the following three algorithms is offered to the user at choice in the considered software :
- Using the classical explicit different scheme;
- Using the classical implicit different scheme;
- Using the method of decomposition of the operator.
The problems of optimum choosing of the algorithm parameters on which depend the accuracy, the time and the possibility of practical realization of the equation solution are realized in these algorithms [5,8].
The method of decomposition of the operator is the original algorithm of solving of the diffusion equation use in multidimensional problems. In this algorithm solution of multidimensional diffusion equation is represented in the form of a linear combination of solutions of some one-dimensional diffusion equations [5,7,8]. Thus, the multidimensional problem is reduced to the one-dimensional one. As the experience shows, the algorithm of decomposition of the operator can be realized by computer much more quickly than the classical explicit and implicit different schemes.
Description of river banks by splines
The problem of the analytical description of plane or spatial region for which the diffusion equations and the boundary conditions are investigated is solved [5,9]. This region represents the part of the channel of the river filled with water in the section for which we are modeling water pollution. The interpolation by splines is used for the analytical description of the plane curve, represented in the form of a sequence of distinct points with given Cartesian coordinates. Such a sequence, in particular, could be one of the river bank lines. The construction algorithms for splines of two types of the simplest explicit form which ensure continuous dependence of the tangent vector of the curve on its parameter are realized in the package. In particular, the river bottom is described by polynomial splines, and the river banks are described by trigonometrical splines [9,10].
Mathematical Models Describing the Pollutants Transport in Rivers
The problem of optimum choice of the steps of digitization of the algorithms of all realized different schemes is solved in the considered software [5,10,11]. The criteria of optimality are the requirements of reduction of the time and the errors of calculation as much as possible. The algorithms are proper for the functions describing the polluting substances transport in water having the derivatives up to the fourth order inclusive.
The total number of nodal points
in the difference schemes determines the time necessary for realization of the algorithm; the accuracy of the obtained result depends on it. Naturally, there arises the question: how to select the numbers
at the given value of
so that the algorithm should be somewhat optimum.
This problem is solved in the following way : at the given values of total nodal points
, there are determined such values of spatial steps of the grid along coordinate axes
for which the upper bound of the module of the residual takes on the minimum value:
is the error of approximation of the equation along the
th coordinate (
The solution of this optimization problem is determined by the formula:
The problem of optimum choice of the value of the parameter
at given spatial steps of the grid is solved in the following way: at the given value of upper bound of the module of the residual of the equation which is equal toε, there are determined such values of
for which the time of calculation takes on the minimum value:
is the error of approximation of the equation by the time coordinate;
is the constant dependent on the algorithm of solution of the sets of linear equations, usually
is the given value of the upper bound of the module of the equation residual.
The solution of this optimization problem is determined by the formula:
The main difficulty for practical realization of the described scheme of optimization is connected with the necessity in estimation of the upper bound of partial derivatives of the function
with respect to the spatial coordinates in terms of which parameters
are expressed without solving the diffusion equation.
Estimation of derivatives of the unknown function
For practical realization of the described optimization schemes, it is necessary to estimate somehow the upper bounds of the modules of partial derivatives of the unknown function with respect to independent variables without solving the initial task. It should be taken into account that the derivatives of these functions can be estimated with the accuracy of the common constant multiplier.
One way of solution of this problem consists in the following: for some values of the parameters
(at solving the diffusion equation), the values of the function
are determined as the first approximation. Then, using the numerical differentiation operations, the required derivatives are determined and the values of the desired parameters are calculated. Using these values, new values of the function
are determined and so on until the difference between the neighboring calculated values of the function are less than the given value .
At solving the diffusion equation, it is possible to use the explicit scheme at the first stage. Then the determination of the values of the sought for function with higher accuracy is necessary. In this case, for calculation of the unknown function values, the iteration method similar to the one used in work  for calculation of the multidimensional integral by the Monte-Carlo method is applied. It is known that the operations of numerical differentiation are not stable, but this should not prevent the realization of the offered method since it requires only rough estimates of unknown derivatives. An obvious drawback of this method is the necessity in performance of a plenty of additional actions. Though, due to the capabilities of modern personal computers, it is negligible at solving the particular tasks.
The effect of smoothness of the inhomogeneous part of the diffusion equation on the accuracy of the results
The inhomogeneous parts of diffusion equation (1) realized in the package contain the impulse functions, which are linear combinations of Dirac delta-functions. These functions and their derivatives are not bounded. Therefore the difference schemes of the solution of diffusion equations are not correct in the vicinities of pollution source localization points.
For elimination of this drawback, it is necessary to consider a model, in which the delta function
is replaced by the bounded numerical function
is the additional parameter. It means that the point sources are replaced by extended ones, the capacity of each of which has the maximum corresponding to the capacity of the source at the point of its location. Such sources can be named quasi-point sources.
should satisfy the following requirements: it reaches the maximum value proportional to
at the point
; it tends to zero at
; its integral between the limits
is equal to unit; its plot represents a bell-shaped curve. The width of such curve is naturally defined as the width of the rectangle which height is equal to the ordinate of the peak of the curve, and its area is equal to the area of the figure limited by the curve and the axis of abscissas. Thus, the parameter
characterizes the width of the plot of the function
, i.e. the sizes of the source. At
, the function
One of the elementary continuous functional dependences which may be used for setting the function
is the trapezoid dependence:
is any real parameter from the interval
. The plot of this function together with the axis of abscissas forms an isosceles trapezium, the top and bottom bases of which have the lengths equal to
, respectively, and the height of which is equal to
. The less is value s, the less this trapezium differs from a rectangle.
Besides the trapezoid dependence, in the developed software package the following dependences are offered to the user’s choice:
In fact, the function
defined by one of the two latter relations plays the part of a smoothing function, which allows us to transform any function from
into an infinitely differentiable function by means of the operation of convolution [10,11].
If the pollution sources do not operate continuously, but they operate for a limited interval of time
, then the capacity of each source should be smoothed not only by spatial coordinates, but also by time. It means that this function should be proportional to
is the additional parameter, and
is the function satisfying the following conditions: at
it tends to unit, and at
to zero. Its plot represents a quasi-stepped curve, which steepness is maximum at
, and this maximum value of the steepness is proportional to
. The parameter
characterizes the time of diminution of the function of discharge. At
, the function
is Heaviside’s stepped function. In the developed software package, it is possible to choose an explicit form of the function
from the following types:
As the computation results show, the smoother are the functions of discharge, the more accurate are the results of solution of the equations in the vicinities of the points of discharge. For given types of functions
, the accuracy of the results increases as the values of parameters