This article will use simplified terminology.
If it is obvious to you without looking at Wikipedia what a matrix, vector and linear approximation of a function are, please understand the fact that most people have long forgotten the university course in higher mathematics, and the article needs to be written in a way that is understandable to everyone.

General calculation principlePreparing data for iterative calculation

Linear approximation of pipeline resistance

Compiling and solving the matrix

For example, this article will use the following scheme:

On the left is the ring network, on the right is its dead-end equivalent, in which the ring is broken and the consumption of consumer 9 is divided into 2 consumers - 27 and 28 (highlighted in cyan). Thus, it is possible to check the correctness of balancing rings in Stokes and other similar programs — the pressure at the inlet to node 5 must match the pressure at nodes 23 and 29.

You can download example stks-file and open it in Stokes. *In the free version of Stokes, select and delete the right half (the equivalent dead-end network), then there will be less than 25 objects and they can be moved and edited. The background with labels can be hidden by pressing* Ctrl+B.

In the case of ring networks (as well as networks with multiple sources), it is impossible to determine in advance how the consumption will be distributed between semi-rings (or between different sources). Therefore, the following sequence of actions is used:

- the flow is distributed in half rings in an arbitrary proportion,
- start of cycle
- the pressure drop is calculated based on these flows,
- the mismatch is calculated at the points of flow merging,
- if the mismatch is within acceptable limits, the calculation is considered completed, exit the cycle
- from overloaded sections, part of the flow is transferred to underloaded sections,
- return to step 2

The mismatch is the difference between the maximum and minimum pressure at the inlet to the point of flow merging from different half-rings/sources. According to SP 42-101-2003, a mismatch of up to 10% is acceptable, but programs for hydraulic calculations (including Stokes) usually continue balancing until it becomes much less than 1%.

The initial flow values are distributed in the same way as in dead-end networks — by sequential passage from each consumer towards the sources. The only difference is that if on the way there is a point of flow merging (which does not happen in dead-end networks), then the flow is divided equally between all the pipes entering it.

In the example network, in the initial state, the point of merging flows is node 4. Its flow of 10 m^{3}/h is distributed equally over 5 m^{3}/h between the pipes from nodes 3 and 5. The initial flows and their directions are shown in orange in the scheme.

In the simplest case, you can simply transfer the entire consumption of an overloaded half-ring to an underloaded one. Then, if necessary, half of this value back, then a quarter in the required direction, etc. until a sufficiently low mismatch is reached. This approach works for networks with one ring and one source.

But for more complex networks, when there are several rings or sources, this approach is not applicable — as long as one ring is balanced, the adjacent one is unbalanced. Therefore, it is necessary to solve the entire system at once. Since the dependence of the pressure drop on the flow through the pipe is nonlinear, the problem is reduced to a system of nonlinear equations, which it is unclear how to solve.

If we reduce the solution to a system of linear equations, then it can be solved, for example, by the Gauss method. Linear equations will look like this:

∆P_{L} = ∆P_{fix} + R Q (1)

где

∆P_{fix} – flow-independent fixed pressure loss, Pa,

R – pipe resistance, Pa / (m^{3}/s),

Q – flow, m^{3}/s.

R and ∆P_{fix} are calculated using the following formulas using the results of previous iterations:

R_{(n)} = (∆P_{L(n-1)} - ∆P_{L(n-2)}) / (Q_{(n-1)} - Q_{(n-2)}) (2.1)

∆P_{fix(n)} = ∆P_{L(n-1)} - R_{(n-1)} Q_{(n-1)} (2.2)

For the first iteration, there are no results from previous iterations yet. Therefore, for the first iteration and for some intermediate ones (more on this below), we calculate R and ∆P_{fix} using the formulas:

R_{(n)} = ∆P_{L(n-1)} / Q_{(n-1)} (3.1)

∆P_{fix(n)} = 0 (3.2)

In these formulas (n), (n-1) and (n-2) denote the iteration number: n are the values used and obtained in the nth iteration. If n=0, we mean the values obtained during the preliminary distribution of flows.

These two options for linear approximation can be represented graphically as follows:

The linear approximation obtained by formulas 2.1 and 2.2 is shown in green; purple – by formulas 3.1 and 3.2; blue – the real dependence of ∆P on Q, calculated by the Darcy-Weisbach equation. Q_{0} – flow after preliminary distribution, Q_{1} – after 1st iteration. The purple R is used for the 1st iteration, the green R and fix are for the 2nd.

For further solution, a square matrix (a system of linear equations) is compiled. The number of equations is equal to the number of pipes in the subnet. Each pipe has a row and a column. The unknowns is Q of the corresponding pipes, the coefficients of the matrix and the vector (column of constant terms) are described below.

For **pipes connected to the consumer**, in the corresponding line all the coefficients are zero except for the coefficient of this pipe - there is 1. In the vector – Q of this consumer: such an entry means that the flow in this pipe is equal to the consumption of the consumer.

In the row **of pipes leading to the splitter** (but not to the point where the flows merge), the coefficients of the corresponding column are 1, the coefficients of the columns corresponding to the pipes flowing from this fork are -1, the rest are 0. Vector 0 is the sum of incoming flows = the sum of outgoing ones.

For **pipes leading to the point of merging of flows**, two types of equations are drawn up:

- For one of the incoming pipes, everything is similar to the usual splitter, only the coefficient of 1 in addition to this pipe is also in the columns of the other incoming pipes.
- For the rest, the equation is compiled using R and ∆P
_{fix(n)}. Two pipe chains are selected from the source to the point of merging of the flows: including the first pipe (from the previous point) and including the current pipe. For each pipe, the coefficient R is set in the first chain, and in the second -R. The vector contains the sum of ∆P_{fix(n)}with opposite signs. The meaning of these equations is that the sum of the pressure losses in the first chain of pipes is equal to the sum in the second.

*In fact, in addition to R and ∆P _{fix(n)}, source pressures and altitude differences are also taken into account, but these nuances are omitted for simplification.*

Further:

- The system of equations is solved by the Gauss method.
- The resulting Q is assigned to the pipes.
- In pipes with Q<0, the flow direction reverses.
- ∆P
_{L}are recalculated using real formulas. - Based on these ∆P
_{L}, R and ∆P_{fix}are calculated for the next iteration.

In addition, the mismatch is calculated at all points of flow merging. If there are no points where it exceeds the permissible threshold, it is considered that the required calculation accuracy has been achieved and the calculation stops. If not, move on to the next iteration. With each iteration, the calculation accuracy increases and the mismatch decreases.