Understanding the research behind Pipe Network Analysis: Hardy Cross Method

Photo by Elimende Inagella on Unsplash

Paper 1

Hardy Cross, developed a mathematical approach for moment distribution in indeterminate structures. He discovered that this method may also be applied for estimating the pressures and discharges in a water distribution network. In the first approach or loop method initial pipe discharges, satisfying the continuity equation at a junction or node, are adjusted to balance the energy equation at a loop. Nodal heads are then obtained from a reference point by adding or subtracting head losses between adjacent nodes.

In the second approach or nodal head method initial heads at nodes are modified in successive iterations to satisfy the continuity equation at a junction. In this method, pipe discharges are estimated by solving the Darcy–Weisbach equations where the head losses are obtained from the head difference between adjacent nodes.

Considering a simple network from H. Cross research paper where h varies with Q².

In the procedure, discharge in certain pipes increase while in others decrease. After repeated correction the answer obtained is inevitably correct because it satisfies the conditions that the quantities balance at each junction, and that the heads balance around each circuit. There are different methods for the hydraulic analysis of water supply networks. In the solution process of most of these methods, a large system of linear equations is solved in each iteration. This usually requires a high computational effort. Hardy Cross method is one of the approaches that do not need such a process and may converge to the solution through scalar divisions. However, this method has two short comings: first, initial discharges should satisfy continuity equation at each node; second a large number of iterations are required to converge to solution. It should be known that the above method was initially developed for hand calculations.

Link for the research paper.

Paper 2

Martin and Peters were the first who used Newton- Raphson method in the analysis of water supply networks[4]. In their method, all equations are written in terms of nodal heads, H. Then the solution is obtained through modifying the head in successive iterations. One of the disadvantages of Martin and Peters’ method is the lack of optimal convergence in large-scale networks. To eliminate this problem, some pipes of the network should be temporarily removed in the analysis procedure. Other disadvantage is the high oscillations to achieve optimal convergence.

Paper: Martin, D. and Peters, G., “The application of newton’s method to network analysis by digital computer”, Journal of the Institute of Water Engineers, Vol. 17, №2, (1963), 115.

Paper 3

In all cases the finite element method was shown to provide superior performance in

comparison to the standard Hardy Cross method. For finite element method to be applied, each

of the element in the pipe network should satisfy the following conditions:

• The algebraic sum of the flows at any joint or node must be zero.
• The value of the piezometric head at a joint or node is the same for all pipes connected to that joint; and
• The flow-head relationship {such as Darcy-Weisbach or Hazen-Williams} must be satisfied for each element or pipe.

The steps for the implementation of the method:

• Taking the input data as number of nodes, elements, known consumption node, known head node, and fluid properties.
• Compute the pipe coefficients, and obtain the initial network characteristics matrix. Then keeping in mind the connectivity matrix, we assemble the matrices into a global matrix.
• Then we apply the boundary conditions and solve.

Advantages of using the above method:

• No need for initial guess.
• Flexibility in applying the boundary conditions.
• Larger and complex networks can be solved.

Link for the research paper.

Paper 4

A pipe network analysis in fluid mechanics analyses fluid flow through a hydraulic network containing multiple or many interconnected branches. The objective is to determine the flow rate and pressure drop of individual network sections. This is a common problem in hydraulic designs. In his paper “Analysis of Streams in Systems of Conduit or Conductor,” published November 13, 1936, Hardy Cross solved the problem of flow distribution in his network of pipes. The inputs and outputs are known, but the flow within the web is unknown. This method can be applied to closed- loop pipe networks (complex sets of parallel pipes). Most engineers commonly use the Hardy Cross method for pipe network analysis. Computer models help engineers solve difficult situations. What is the maximum fire current at a given point in the system? How long can this flow of fire be provided? Do you need pipes? Will there be adequate pressure and flow if a residential or commercial building is built? If not, what lengths and sizes of water pipes will developers need to upgrade to accommodate the proposed construction?

Link to the research paper.

Paper 5

Piping networks have many practical applications, from water and gas distribution systems to air conditioning. Simple problems, such as a single branch connecting two reservoirs, can be solved analytically, but complex network problems require an iterative approach using digital computers. The Hardy-Cross method is the most common method for solving this type of problem. The Hardy Cross method is an iterative method for determining flows in pipe network systems where the inlet and outlet flows are known, but the flows within the network are unknown. The pipe’s length, diameter, roughness, and other essential characteristics must be known. Before the introduction of this process, complex distribution piping systems were complicated to solve due to the non-linear relationship between pressure drop and flow rate. The water distribution system model has become an essential and valuable tool for civil engineers. Models are often used to optimize the design of new distribution systems or to analyze significant expansions or modifications to existing distribution systems. The introduction of the Hardy Cross method of pipe flow network analysis revolutionized municipal water supply design.

• What is the maximum fire flow at a given point in the system?• How long can that fire flow be provided for?
• What size pipe installation would be necessary between two points in a system to increase the pressure at one of the points to the minimum force required?
• If a subdivision or commercial development is built, will adequate pressures and flows exist? If not, what length and size of water mains must be upgraded by the developer to allow for the proposed construction?

Link for the research paper.

Problem Statement

Consider a pipe network comprised of nine (J = 9) junctions, thirteen (N = 13) branches and four (M = 4) loops. The data provided include: the external inflows (Qin) and outflows (Qout) at all network junctions, and the resistance coefficients (Ri) and pipe diameters (Di) for all branches. For this given network, write a computer code for the determination of the discharge through each branch pipe following Hardy-Cross method. Set the convergence criterion for the maximum correction value of ΔQl as 0.001.

1. main.m

clc
pipe_network = load('pipe network.txt');
connectivity_matrix = load('connectivity matrix.txt');
branch_resistance = load('branch resistance.txt');
br_res = branch_resistance(:,3);
junctions = pipe_network(1,1);
loops = pipe_network(2,1);
branches = pipe_network(3,1);
Q_in = pipe_network(4,1);
Q_out = Q_in - (8*55);

% There are four loops in the question
% Asumming the data in the pipelines following the continunity principle
%disp(Q_out);
br_flow = (assume_data(Q_in))';

% Finding the correction
br_flow_corrected = flow_correction(br_res,br_flow);
data = table(br_flow,br_flow_corrected);
disp(data);

2. flow_correction.m

function [br_flow] = flow_correction(br_res,br_flow)
    sigma_1 = [-1,1,0,0,1,0,0,1,0,1,0,0,0];
    sigma_2 = [0,0,1,1,-1,-1,0,0,0,0,0,0,0];
    sigma_3 = [0,0,0,0,0,1,1,-1,-1,0,0,0,0];
    sigma_4 = [0,0,0,0,0,0,0,0,1,-1,1,1,1];
    loop1 = [1,2,5,8,10];
    loop2 = [3,4,5,6];
    loop3 = [6,7,8,9];
    loop4 = [9,10,11,12,13];
    
    while true
        [del_q1,del_q2,del_q3,del_q4,br_flow] = del_q(br_res,br_flow,sigma_1,sigma_2,sigma_3,sigma_4,loop1,loop2,loop3,loop4);
        if del_q1 < 0.001 && del_q2 < 0.001 && del_q3 < 0.001 && del_q4 < 0.001
            break;
        end
    end
end

3. del_q.m

function [del_q1,del_q2,del_q3,del_q4,br_flow] = del_q(br_res,br_flow,sigma_1,sigma_2,sigma_3,sigma_4,loop1,loop2,loop3,loop4)

         num = 0;
         den = 0;
         for i = loop1
            num = num + sigma_1(1,i) * br_res(i,1) * br_flow(i,1) * abs(br_flow(i,1));
            den = den + 2 * br_res(i,1) * abs(br_flow(i,1));
        end
        del_q1 = num/den;
       for i = loop1
           
           if abs(del_q1) < 0.001
               break;
           end
           
           br_flow(i,1) = br_flow(i,1) - sigma_1(1,i) * del_q1;
      end
        num = 0;
        den = 0;
       for i = loop2
           num = sigma_2(1,i) * br_res(i,1) * br_flow(i,1) * abs(br_flow(i,1)) + num;
           den = 2 * br_res(i,1) * abs(br_flow(i,1)) + den;
       end
       del_q2 = num/den;
       for i = loop2
           if abs(del_q2) < 0.001
                break;
           end
           br_flow(i,1) = br_flow(i,1) - sigma_2(1,i) * del_q2;
       end
        num = 0;
        den = 0;
       for i = loop3
           num = sigma_3(1,i) * br_res(i,1) * br_flow(i,1) * abs(br_flow(i,1)) + num;
            den = 2 * br_res(i,1) * abs(br_flow(i,1)) + den;
        end
        del_q3 = num/den;
        for i = loop3
           if abs(del_q3) < 0.001
                break;
            end
            br_flow(i,1) = br_flow(i,1) - sigma_3(1,i) * del_q3;
        end

        num = 0;
        den = 0;
       for i = loop4
            num = sigma_4(1,i) * br_res(i,1) * br_flow(i,1) * abs(br_flow(i,1)) + num;
            den = 2 * br_res(i,1) * abs(br_flow(i,1)) + den;
        end
        del_q4 = num/den;
        for i = loop4
            if abs(del_q4) < 0.001
                break;
            end
            br_flow(i,1) = br_flow(i,1) - sigma_4(1,i) * del_q4;
        end
end

4. assume_data.m

function [br_flow] = assume_data(Q_in)
    % We assume the flow in clockwise direction with positive sign.
    br_1 = Q_in * rand(1,1);
    br_2 = Q_in - br_1;
    br_5 = (br_2 - 55) * rand(1,1);
    br_3 = br_2 - 55 - br_5;
    br_4 = br_3 - 55;
    br_8 = (br_5 - 55) * rand(1,1);
    br_6 = br_5 - 55 - br_8;
    br_7 = br_6 + br_4 - 55;
    br_9 = (br_8 - 55) * rand(1,1);
    br_10 = (br_8 - 55) - br_9;
    br_13 = br_7 + br_9 - 55;
    br_12 = br_13 - 55;
    br_11 = br_12 - 55;
    br_flow = [br_1,br_2,br_3,br_4,br_5,br_6,br_7,br_8,br_9,br_10,br_11,br_12,br_13];
end