Mathematical modeling of fluid flow in complex multi-channel structures

The work is focused on a numerical hydraulic model, which allows one to simulate fluid flow in a complex, arbitrarily configured multi-channel system with a capability of automatic visualization of its structure. The main elements of a system are channels which are united in the system by means of a particular set of coupling elements (local resistance, knots). The processes models in channels are based on the equation of continuity, momentum conservation. The system structure is made by forming matrices of regional and internal boundary conditions according to the developed algorithm. By the model the algorithm is created, on the basis of which the application program is developed. The application program enables one to determine parameters of a steady state of fluid flow in the complex multichannel structures. The calculation structure is shown.


Introduction
In the analysis and model operation of the complex hydraulic systems the development of methods, models, and analysis algorithms of the fluid flow in channel systems with an automated display of the structure of systems is urgent [1][2][3][4]. Realization of such an approach in the form of the applied program allows the creation of a universal tool for the analysis of fluid flow in the channel systems of any given structure [5][6][7][8][9].
For calculation of the fluid flow in the complex multi-channel systems, an approach for the calculation of the steady state of flow in the multi-channel system of any given structure with an automated display of the structure of system is realized. It relies on the representation of the mathematical model of a system in the form of the homogeneous boundary value problem and the automated formation of matrices of boundary conditions. This approach is more simple in comparison to the existing methods of model operation and the analysis of processes in the complex hydraulic systems, but at the same time allows investigating them efficiently. Relying on this approach, the computing model of calculation of the steady state of flow in the complex multi-channel systems of any given structure is constructed [10].

The mathematical model
It is given that channels, local resistance, and connecting clusters of channels can be elements of the system.
The system of one-dimensional equations of continuity, momentum conservation describing fluid flow in the channel is presented below: ṁ is mass flow of an actuating fluid in the channel, kg/s; P is pressure of an actuating fluid in the channel, Pa; i is enthalpy of an actuating fluid in the channel, J/kg; ρ is density of an actuating fluid in the channel, kg/m 3 ; F is a sectional area of the channel, m 2 ; П w is a perimeter of transverse section of the channel, m; τ w is a tangential stress of friction on a wall, N/m 2 . The set of equations (1) -(2) is supplemented with the closing dependences from friction: ξ is a coefficient of hydraulic resistance.
For channel connection knots equations corresponding to a mass conservation law (4), and conditions of pressure equality in the knot are applied (5): indexes: j is a channel of supply, i is a diverting flume. For local resistances, it is considered that their hydro-resistance depends on expense, pressure in initial sections. It is offered to consider that the local resistance similar to a throttling orifice is a simple knot which connects 1st (j) supplying-and 2nd (i) allocating-channels. The equations which describe the processes on the throttling orifice have the following forms: The description of interaction of the system with a surrounding medium requires boundary conditions at the system input and output, i.e. for those points where fluid comes into the system and leaves it (supply and output devices).
For entry section of a supplying channel (index n) conditions are accepted: For an allocating channel (index k): The ratios (7) and (8) must be written down for all heating mains (supplying and leaving), which can be a few in general case.
Hydraulic characteristics of local and induced resistance are: K is a loss coefficient. The shown equations and ratios written down for each element of the system form a self-contained model of the system's steady state. The model of steady conditions of the system represents a non-linear boundary value problem. For its solution, it is necessary to construct a convergent repetitive computing process.
Applying the procedure of quasi-linearization to the equations (1) -(2) and supplementing them with relaxation parameters, the following equations are received: n is the number of iteration; t is a relaxation parameter.
For the automated development of the model of steady-state conditions the state vector Y is added: ,..., N is the number of channels in the system Then the set of equations (10) -(11) can be written down in the following matrix form: where the matrix Θ is 2N×2N in size and the vector Φ is constructed of coefficients of the equations (10)-(11). Thus, Θ and Φ include parameters only from the iteration n, therefore, the set of equations (13) is linear with respect to the components of the vector Y.
The linearization of boundary conditions gives the following ratios: −for the channel connection knots: − for the input supply channel: − for the output tail-race channel: After the linearization a set of boundary conditions for the circulating multi-channel system consisting of the elements stated above can be presented in the following matrix form: where B and C are matrices with 2N×2N in size, and D is the column vector containing 2N elements. The matrices B, C, D are formed on the basis of equations (14) − (17), taking into account specified system structure, the stationary state of which is required to be defined. The solution to the boundary value problem (13), (18) until reaching the repetitive process convergence will give the required set of parameters (pressure, mass flow, an enthalpy) of the steadystate of the system. This solution is obtained based on a sweep method (reduction to a Cauchy problem).
The matrices B, C, D contain full information on the structure of the system. Thus, the automated design of the computing steady-state model is reduced to the formation of the specified matrices by a special algorithm.
The algorithm of formation of matrices B, C, D consists in the following. Each of the ratios (14) − (17) can be presented in the matrix form: where j is the number of system knot. Each of the matrices , has 2N columns and M lines. The vector has M lines, where M is the number of boundary conditions generated by this element. Formation of these matrices is made in accordance with the structure of the equations (14) -(17).
The scheme below shows that the filling of the matrices B, C, D from (18) is made on the block principle: where k is the number of knots in the system.
To implement the sweep method, the problem (13), (18) is reduced to the problem: For this purpose, the vector is supplemented with one more component y 2N+1 , the system of differential equations (13) is supplemented with the equation the matrix A is supplemented with the zero line, the column Φ and the element a 2N+1, 2N+1 = 0, and boundary conditions are supplemented with the condition y 2N+1 (x 2N+1 = 0) = 1. The sweep method (reduction to a Cauchy problem) consists in the following: Let us suppose it is necessary to find the solution of combined equations which meets the boundary conditions Then in the point x=0 If matrix C(x) is known, then, by solving system (24), it is possible to determine initial condition Y(0) for combined equations (21).
Differentiating (23) and considering (21), the following equation is received: As the last equality must be fulfilled for Y(x) at any value of x from the interval [0, L], then where A T and C T are a transpose of matrices A and C. Therefore, in order to determine matrix C(0)=(C T (0)) T , it is necessary to integrate combined equations (25) from x=L to x=0 with the initial condition of C T (x)≡C T , as condition (23)  System (25) can be written down as follows: where is the i-th column of matrix C T (x) coinciding with the i-th row of matrix C(x). Besides, as system (29) shows, if the k-th row of matrix C(x)≡C consists of zeroes, then matrix C(x) has the k-th row consisting of zeroes throughout the interval [0, L].
After matrix C(0) is found, the system (24) for finding initial condition Y(0) is used. Then the Cauchy problem is integrated for combined equations (21), with the determined initial condition from the point 0 = x to the point x=L. At the same time, the solution for the boundary value problem (21), (22) is found.

Algorithm and application program
Based on the model described above, an algorithm, serving as a foundation to create an application program which allows determining parameters (the size of mass flow, pressure in channels) of the stationary state of a system of channels in any given structure, is developed [ The structure of the studied system is described. Channels and knots are numbered sequentially from the entrance into the system to its exit. For each channel, its diameter, length and friction resistance coefficient are determined. Characteristic of each knot (local resistance, entrance into the channel, exit out of the channel, prime knot of connection of channels), the numbers of channels for which this knot is the beginning, and the numbers of channels for which this knot is the end is described. Depending on the knot's characteristic, the pressure of medium at the entrance of each of the approach channels and the outlet pressure of each of the diversion channels is set.
The program forms matrices of boundary conditions. On the basis of each knot's description, the necessary elements of these matrices are filled. It is carried out by a consequential, following the numbering of knots, calculations of the filled lines of matrices, depending on the equations that describe a specific knot and on the numbers of channels that are included in this knot of calculations of the filling columns.
The sweep method numerically solves the homogeneous boundary value problem. Computational solution of the system of ordinary differential equations of first order with given initial conditions (the Cauchy problem) is carried out by the Runge-Kutta method.
The values of mass flows in the system's channels, obtained during the current iteration, are compared to mass flows from the previous iteration. If divergence is more than the given accuracy, then matrices of boundary conditions are formed again with the obtained values, and the calculation is repeated. The iteration process is continued till convergence.

Conclusion
The simulation results allow for quick and visual investigation of complex multi-channel structures' operation at a design stage and during its exploitation. The scheme of channels' connection can be of any complexity.