Fluid Mechanics & Heat Transfer

Textbook Solutions: Fluid Mechanics, Heat & Mass Transfer, Aerodynamics

FM, HTX, MTX, AERO, COMBUSTION

This page is being continuously updated with complex (text-book type) problems in Fluid Mechanics, Heat Transfer, Aerodynamics, Mass Transfer, Combustion and Thermodynamics. The solutions to compressible flows including sub-sonic, sonic and supersonic flows inside a converging-diverging nozzle will be presented.

Table of Contents:

Summary of Hydrostatics

• For a plane surface immersed in a liquid in any orientation, the total normal pressure over the surface is P = A * ρ * g * h where h is the vertical depth of the centre of the gravity of the surface from free surface of the liquid, A is area of the surface.
• A plane surface immersed in a liquid is acted upon by a system of parallel forces which can be reduced to act on a single point of the surface known as centre of pressure.
• A curved surface immersed in a liquid is not acted upon a "system of parallel forces" and hence in general they cannot be reduced to a single force.

Centre of pressure of vertical rectangular plate

This is a very useful case which has frequent engineering applications.

Similarly, for a triangular plate with base horizontal and at depth of Y₂ and vertex up at depth of Y₁, the centre of pressure is given by:

Centre of pressure of vertical circular plate

Even though the complicated definite integration is possible analytically, this is the time the introduce a simple method based on properties of sections. Note that the numerator of expression for centre of pressure is the "Second Moment of Inertia" or simply the Moment of Inertia about a "chosen axis" - typically the free surface of the liquid. The denominator is the "Moment of Area" about the same axis as chosen above - also known as "Statical Moment". Thus:

• If the area has axis of symmetry about the vertical axis, the centre of pressure will fall on it.
• Conversely, the centre of pressure will fall on the medial line, the vertical axis passing through the centre of gravity on the horizontal axis.
• The position of the centre of pressure is always below the centre of gravity of a surface.

Buoyancy - Example Problem

A metallic block is floating along with a wooden block as shown in the figure below. If this metallic block falls into the tank, calculate the new depth of floatation and height of water level in the tank.

Analogy Between Fluid Flow, HTX and MTX

The method of relating two phenomena with some come feature is a powerful tool to both understand and remember the concepts. Following table summarizes key parameters which are analogous in the field of fluid flow (momentum transfer), heat transfer and mass transfer.

1D Heat Transfer with Convective Boundaries and Heat Generation

The designs of heat exchangers are the most predominant application of conductive heat transfer with such boundary conditions. The temperature profile in walls of a water-cooled internal combustion engine is also subject to this combination of boundary conditions though heat generation is not present in this application. Some applications can be current carrying conductor cooled by convection on both inner and outer diameters, nuclear fission in a annular cross-section cooled by convection on both the inner and outer radii. The governing equation and general solution of the differential equation is given by

Heat flow is assumed positive from left to right. Hence, the convective heat transfer on left face is negative if ambient temperature is < wall temperature on the left face!

Plug Flow between Parallel Plates

Flow in an Annular Passage

Flow through Pipe Branches forming Tee or Wye Network

Note velocities in pipe branches 1-2, 2-3 and 3-4 are not known in advance and hence the loss coefficients K_MN cannot be adjusted to same velocity. Note that the junction losses at node 2 can be incorporated either in K₁₂ or in K₂₃ and K₂₄. It is more convenient to club the node loss in K₁₂ as appropriate assignment into K₂₃ and K₂₄ will involve extra calculations.

This is a non-linear equation in V₂₃ and needs to be solved using trial-and-error or iterative approaches. Microsoft Excel Goal Seek utility can be used to solve this equation as well. A good initial guess would be required.

A more compact approach in terms of fluid resistances can be used as demonstrated below.

Thus, we have 3 equations and 3 unknowns – but they are still non-linear and needs to be solved using iterative method.

Force on a curved pipe with varying cross-section

A reducing bend having diameters D upstream and d downstream, has water flowing through it at the rate of m [kg/s] under a pressure of p1 bar. Neglecting any loss is head for friction calculate, the horizontal and vertical components of force exerted by the water on the bend.

Newton's Law of Motion applied to Fluid Elements:

Newton's Third Law Applied to Pipe:

This is the force acting on the fluid and only means to apply the force on the fluid is the walls of the pipe. Hence, from Newton’s third law, and equation and opposite force will act on the pipe.

Force on fluid element due to pressure at boundaries

If pipe is straight

Equation of State: Soave-Redlick-Kwong Equation

or gases at very high pressure and those in liquid state such as Liquified Petrolium Gas (LPG), the accuracy of ideal gas equation falls sharply. Many improvizations have been made and SRK equation is one such method. The approach is demonstrated using following sample problem:

A gas cylinder with a volume of 5.0 m³ contains 44 kg of carbon dioxide at T = 273 K. estimate the gas pressure in [atm] using the SRK equation of state. Critical properties of CO₂ and Pitzer acentric factor are: Tc = 304.2 K, pc = 72.9 atm, and ω = 0.225.

General cubic equation of state is given by following equation

Alternatively, the equation can be solved using Newton's method as described below.

dp = ρ . {F_X dx + F_Y dy + F_Z dz + ω²(xdx + ydy)} where F_X is the force per unit volume is X-direction. Alternatively, dp = ρ . {F_X dx + F_Y dy + F_Z dz + ω²rdr}. Any point where liquid separates is the location where pressure vanishes i.e. p = 0. Latus rectum of free surface = 2g/ω².

Half Filled Rotating Tube: A circular tube is half full of liquid and is made to revolve around a vertical tangent line with angular velocity ω. If radius of the tube is R, calculate the angle between horizontal and diameter passing through the free surface of the liquid.

The tube is rotating about red line and the liquid level is shown by yellow zone. Point 'D' is the lower most point of free surface, θ is the angle which needs to be calculated.

The pressure at any location is given by: dp = ρ (ω²r.dr + gdz) where r is the distance of the point from the axis of rotation and z is the distance of the point from the origin. The point 'D' is rotating at radius r = R(1-cosθ), z = R.sinθ and is at p = 0. Another known point 'C' is at pressure zero, r = R(1+cosθ), z = -Rsinθ. Note that +z is along the downward direction (this is only a convention).

Liquid filled in rotating conical pot: A conical pot of height H and vertical (full cone) angle θ contains liquid k-times (x < 1) the volume of full cone. Find out the minimum angular velocity at which shall prevent liquid overflowing the cone.

Thus; z = (1-k).2H/3 ... [1] Also, from paraboloid of revolution (latus rectum): 2g/ω² = R²/z ⇒ z = ω².R²/2/g ... [2]

From equations [1] and [2], ω².R²/2/g ≤ (1-k).2H/3 ⇒ ω ≤ [4g/3.(1-k).H/R²]^1/2

Fluid Pressure in Rotating Reference Frame

For a container rotating cylindrical about its axis: the shape of the free surface is a parabola and fluid inside the rotating cylinder forms a paraboloid of revolution, whose volume is one-half of the volume of the "circumscribing cylinder". To calculate angular velocity at which the liquid at the center reaches the bottom of the cylinder just as the liquid at the curved wall reaches the top of the cylinder: ω_spill = (2gH)^0.5/R.

Ball in rotating tube

The centroid of a hemi-sphere is at 3R/8 from the base. Using this value, pressure on the left half of the ball = 1/2.ρ_F ω²(r-3R/8)². The pressure on the right half of the ball = 1/2.ρ_F ω²(r+3R/8)². Net force acting on the ball towards the axis of rotation = 1/2.ρ_F ω²(4.r.3R/8)² × π*R² where projected area = π*R².

Net force due to fluid pressure towards axis of rotation, F_F = 1/2.πR².ρ_F ω²(r.3R/2) = 3/4πR³.ρ_F ω²r

Centrifugal force due to own mass of the ball, F_B = 4/3.πR³.ρ_B ω²r

A circular cylinder of radius 'a' whose axis is vertical, is filled to depth h with homogeneous liquid of density ρ. Find the pressure as function of radial and vertical positions. Also, find the forces acting on the top lid if closed.

dp = ρ{-g.dz + ω²rdr} -> p = ρ{-g.z + 1/2.ω²r²} + C

A hollow cone nearly filled with water, rotates uniformly about its axis which is vertical, the vertex being uppermost. Find the total thrust on the base of the cone. 'a' is the radius of the base and 2α the vertical angle of the cone.

p = ρ{-g.z + 1/2.ω²r²} + ρgh

Practice Questions:

1. A hollow sphere of radius 'R' is half-filled with liquid and is rotating at angular velocity ω about its vertical diameter. Find out the value of ω at which the lowest point of the sphere is just exposed. Answer: [2g/R]^1/2.[2-4^1/3]^-1/2. Hint: volume of water remains constant.
2. A narrow annular horizontal tube of radius R containing liquid is filled to half of its circumference and is rotating about the vertical diameter with angular velocity ω. Calculate the value of angular velocity at which liquid will just separate. Answer: [2g/R]^1/2. Hint: pressure will become zero at lowest point when liquid separates.
3. A vertical inverted conical vessel (base at the top) of height 'h' and cone angle 2θ is half-filled with a liquid oil. Find out the maximum angular velocity ω at which the oil shall not overflow. Answer: [2g/3h]^1/2.cot(θ) Hint: use latus rectum as defined earlier.

Steps to Create a Thermal Network Solver

Step-0: Define problem set-up: Steady or Transient, Initial Temperature, Total Run Time, Time Step Solution, Time Step Data Save.

Step-1: Define input table for the network. For Steady state cases, Cp and Mass shall be neglected.

Note that: NND = number of nodes, NCV = Number of control volumes = NND+1, XND[0] = XCV[0] = 0, XCV[i] = (XND[i+1] + XND[i] ) / 2 for i = 1 to NND-1, XCV[NCV] = XND[NND]. LCV[0] = XND[1]/2, LCV[i] = XCV[i] - XCV[i-1], LCV[NCV] = 0.5 x (XND[NND] - XND[NND-1])

Step-2: Define resistance table. All information can be supplied as CSV file in following order (use negative number for boundary nodes). Here, Tref can be used to specified reference temperature for nodes with HTC boundary conditions as well as nodes with known or fixed temperatures. Value of 0 for Temperature or Tref indicate unknown or optional. There are some limitations to define the nodes such as a node must be placed where there is step change of cross-section. The node numbering must start with 1 and increase sequentially without any gap.

           [3]      [6]
            |        |
            |        |
[1]--------[2]------[5]-------[8]-------[9]-------[10]
            |        |
            |        |
           [4]      [7]
R_ID, s_node, e_node, xLen,  xAr,      k,     Cp,    Den 
R01,  1,      2,      0.050, 0.0005, 200.0, 500.0, 2700.0
R02,  2,      3,      0.025, 0.0005, 200.0, 500.0, 2700.0
R03,  2,      4,      0.025, 0.0005, 200.0, 500.0, 2700.0
R04,  2,      5,      0.125, 0.0005, 200.0, 500.0, 2700.0
R05,  5,      6,      0.025, 0.0005, 200.0, 500.0, 2700.0
R06,  5,      7,      0.025, 0.0005, 200.0, 500.0, 2700.0
R07,  5,      8,      0.025, 0.0005, 200.0, 500.0, 2700.0
R08,  8,      9,      0.025, 0.0005, 200.0, 500.0, 2700.0
R09,  9,     10,      0.025, 0.0005, 200.0, 500.0, 2700.0

def calcThermalRes(xL, xA, xk):
  return xL/xA/xk
def calcThermalIntertia(xL, xA, xDen, xCp):
  return xDen * (xL * xA) * xCp

NN stands for Node Number
-------------------------
NN,   Qdot,  HTC,    Tref
 1,    0,     0,     288
 2,   50,     0,     0
 3,   50,     0,     0
 4,   50,     0,     0
 5,    0,     0,     0
 6,    0,     0,     323
 7,    0,    10,     298
 8,    0,    10,     298
 9,    0,     0,     0
10,    0,     0,     373
 

import csv
def csvToColumnLists(csv_file):
  '''
  Reads a CSV file and returns a dictionary where keys are column headers
  and values are lists containing the data for each column. First non empty
  row is assumed to be header containing strings without spaces. Header names
  are case-sensitive.
  '''
  column_data = {}
  with open(csv_file, 'r', newline='') as file_csv:
    reader = csv.reader(file_csv)
    headers = next(reader)
    headers = [header.strip() for header in headers]
    # Initialize empty list for each header and append data
    for header in headers:
      column_data[header] = []
    for row in reader:
      for i, value in enumerate(row):
        column_data[headers[i]].append(value.strip())
  return column_data
#
data_r = csvToColumnLists('thermal_r.csv')
node_col_1 = data_r['s_node']
node_col_2 = data_r['e_node']
node_col_1 = [int(NN) for NN in node_col_1]
node_col_2 = [int(NN) for NN in node_col_2]
xL = data_r['xLen']

# Find names of nodes as union of sets from column 2 and 3
NNU = list(set(node_col_1) | (set(node_col_2)))
NNU = [int(N) for N in NNU]
NNU.sort()

NNL = node_col_1 + node_col_2
max_n = int(max(set(NNL), key=NNL.count))
NNR = len(data_r['R_ID'])
NMAT = [[0] * max_n] * NNR

Create thermal resistance 2D array

def createRTH(xLen, xAr, xk, node_i, node_j):
  '''
  Create thermal resistance as 2D matrix, the size of matrix is 1 greater than
  the maximum node number in the network. Note that node numbering starts from
  1 and increase sequentially without any gap or step.
  '''
  NTH = int(max(max(node_i), max(node_j)))
  NR = len(node_i)
  RTH = []
  for _ in range(NTH):
    row = [0] * NTH
    RTH.append(row)
    
  for i in range(0, NR):
    m = int(node_i[i]) - 1
    n = int(node_j[i]) - 1
    RTH[m][n] = float(xLen[i]) / float(xAr[i]) / float(xk[i])
    RTH[n][m] = RTH[m][n]
  return RTH

Step-3: Define governing equations using energy balance at each nodes with unknown temperatures and nodes with known temperatures.

For node-1 above: (T₁ - T₂)/R₁₂ + (T₁ - T₃)/R₁₃ = Qdot₁

or

A₁ T₁ + B₁₂T₂ + B₁₃T₃ = Qdot₁ where the second values in the subscripts correspond to entry under "Connected Nodes".

For node-2 above: (T₂ - T₁)/R₁₂ + (T₂ - T₃)/R₂₃ + (T₂ - T₄)/R₂₄ = Qdot₂

or

A₂ T₂ + B₂₁T₁ + B₂₃T₃ + B₂₄T₄ = Qdot₂ where the second values in the subscripts correpond to entry under "Connected Nodes".

Step-4: Convert the individual equations into equivalent index notation (subscripts converted into indices of arrays).

Step-5: Define variables and arrays. NND = number of nodes, 1D arrays: XND[NND], XCV[NCV] Qdot[NND], Cp[NND], mass[NND], cond[NND], XAR[NND], LCV[NND]. 2D arrays: coeff_a[NND, NND], coeff_b[NND, NND], coeff_c[NND, NND], th_res[NND, 10], con_nodes(NND, 10). Each node is assumed to be connected to maximum 10 nodes.

Write functions to approximate differential equations

For steady state 1D heat conduction with volumetric heat generation and convective heat loss:

Another example from web:

Training Project: Development of thermal network model and validation of results with 3D CFD simulations using ANSYS Fluent or OpenFOAM is a good internship project. This will impart not only programming skills but also help the engineer gain insight into heat transfer and matrix algebra.

Thermal Network Example Hand Calculations

q'[100W]----->[1]-------[2]--------[3]== T [300 K]
                  L: 0.10 [m]
                  A: 5E-4 [m2]
                  k: 50.0 [W/m-K]

Node-1: 
 1/R12 * [T2 - T1]                     = q'
Node-2:
 1/R12 * [T2 - T1] + 1/R23 * [T2 - T3] = 0
-1/R12 * [T1] + [1/R12 + 1/R23] * [T2] = T3/R23
-1/R12 * [T1] + [1/R12 + 1/R23] * [T2] = T3/R23

Matrix to solve

-1/R12 * [T1] + 1/R12           * [T2] = q'
-1/R12 * [T1] + [1/R12 + 1/R23] * [T2] = T3/R23

Some other configurations that can serve as validation for any thermal network solver are:

q'[100 W]----->[1]-------[2]--------[3]== T [300 K]
                          |
                          |
                         [4]
                      q' [50 W]


                      HTC: 20 [W/m2-K], TREF=300 [K]
                         [5]
                          |
                          |
q'[100 W]----->[1]-------[2]--------[3]== T [300 K]
                          |
                          |
                         [4]
                       q' [50 W]

Generic node in a thermal network:

                [s]   [r]   [p]
              .'  \    |   /
             '     \   |  /
            '        *[n]'--------  [m]
            '          |
             ` ,       |
                 ` - =[u]

[Tn - Tm]/RTH[n,m] + [Tn - Tp]/RTH[n,p] + ... + [Tn - Tu]/RTH[n,u] + h.A.[Tn - Tref] = q

Σ_j(1/RTH[n, j]) . Tn - {Σ_j[1/RTH[n,j] . Tj} = q + h.A.Tref ... ... ...[1]

Σ_j(CTH[n, j]) . Tn - {Σ_j[CTH[n,j] . Tj} = q + h.A.Tref ... ... ...[2]

Next step: create thermal resistance 2d array. RTH[j,k] are calculated earlier and known for all j and k - sparse symmetric matrix, RTH[j, k] = RTH[k, j] and diagonal entries = 0. Note that thermal resistances by definition are always positive, a negative value is used to indicate direction of heat transfer.

Thermal Resistance Array:
       1 |     2 |     3 |     4 |     5 |     6 |     7 |     8 |     9 |    10 |
 --------|-------|-------|-------|-------|-------|-------|-------|-------|-------|
 1     0 |  0.50 |     0 |     0 |     0 |     0 |     0 |     0 |     0 |     0 | 
 2  0.50 |     0 |  0.25 |  0.25 |  1.25 |     0 |     0 |     0 |     0 |     0 | 
 3     0 |  0.25 |     0 |     0 |     0 |     0 |     0 |     0 |     0 |     0 | 
 4     0 |  0.25 |     0 |     0 |     0 |     0 |     0 |     0 |     0 |     0 | 
 5     0 |  1.25 |     0 |     0 |     0 |  0.25 |  0.25 |  0.25 |     0 |     0 | 
 6     0 |     0 |     0 |     0 |  0.25 |     0 |     0 |     0 |     0 |     0 | 
 7     0 |     0 |     0 |     0 |  0.25 |     0 |     0 |     0 |     0 |     0 | 
 8     0 |     0 |     0 |     0 |  0.25 |     0 |     0 |     0 |  0.25 |     0 | 
 9     0 |     0 |     0 |     0 |     0 |     0 |     0 |  0.25 |     0 |  0.25 | 
10     0 |     0 |     0 |     0 |     0 |     0 |     0 |     0 |  0.25 |     0 |

       1 |     2 |     3 |     4 |     5 |     6 |     7 |     8 |     9 |    10 |
 --------|-------|-------|-------|-------|-------|-------|-------|-------|-------|
 1     0 |   2.0 |     0 |     0 |     0 |     0 |     0 |     0 |     0 |     0 | 
 2   2.0 |     0 |   4.0 |   4.0 |   0.8 |     0 |     0 |     0 |     0 |     0 | 
 3     0 |   4.0 |     0 |     0 |     0 |     0 |     0 |     0 |     0 |     0 | 
 4     0 |   4.0 |     0 |     0 |     0 |     0 |     0 |     0 |     0 |     0 | 
 5     0 |   0.8 |     0 |     0 |     0 |   4.0 |   4.0 |   4.0 |     0 |     0 | 
 6     0 |     0 |     0 |     0 |   4.0 |     0 |     0 |     0 |     0 |     0 | 
 7     0 |     0 |     0 |     0 |   4.0 |     0 |     0 |     0 |     0 |     0 | 
 8     0 |     0 |     0 |     0 |   4.0 |     0 |     0 |     0 |   4.0 |     0 | 
 9     0 |     0 |     0 |     0 |     0 |     0 |     0 |   4.0 |     0 |   4.0 | 
10     0 |     0 |     0 |     0 |     0 |     0 |     0 |     0 |   4.0 |     0 |

Next step: find out nodes with known temperature.

data_n = csvToColumnLists('thermal_n.csv')
node_nnum = data_n['NN']
node_nnum = [int(NN) for NN in node_nnum]
node_hgen = data_n['Qdot']
node_hgen = [float(x) for x in node_hgen]
node_htc  = data_n['HTC']
node_htc  = [float(x) for x in node_htc]
node_Tref = data_n['Tref']
node_Tref = [float(x) for x in node_Tref]
node_known_T = [
  item_1 for item_1, item_2, item_3 in zip(node_nnum, node_htc, node_Tref)
  if item_2 == 0 and item_3 != 0
]

Node connectivity array:
    2 |     1 |     3 |     4 |     5 | 
    3 |     2 |     0 |     0 |     0 | 
    4 |     2 |     0 |     0 |     0 | 
    5 |     2 |     6 |     7 |     8 | 
    7 |     5 |     0 |     0 |     0 | 
    8 |     5 |     9 |     0 |     0 | 
    9 |     8 |    10 |     0 |     0 | 

           [3]      [6]
            |        |
            |        |
[1]--------[2]------[5]-------[8]-------[9]-------[10]
            |        |
            |        |
           [4]      [7]

Using Σ_j(CTH[n, j]) . Tn - {Σ_j[CTH[n,j] . Tj} = q + h.A.Tref ... ... ...[2]

(C[2, 1] + C[2, 3] + C[2, 4] + C[2, 5]) * T[2] 
                              - C[2, 3] * T[3] 
                              - C[2, 4] * T[4] 
                              - C[2, 5] * T[5] 
                              = C[2, 1] * T[1] + q[2] + h[2] * A[2] * (T[2] - Tref[2])

A[0, 1]*T[2] + A[0, 2]*T[3] + A[0, 3]*T[4] + A[0, 4]*T[5] = B[1]

C[3, 2] * (T[3] - T[2]) = q[3] + h[3] * A[3] * (T[3] - Tref[3])

A[1, 0]*T[2] + A[1, 2]*T[3]                               = B[2]

Node-4: convection node

C[4, 2] * (T[4] - T[2]) = q[4] + h[4] * A[4] * (T[4] - Tref[4])

A[2, 0]*T[2]                + A[2, 3]*T[4]                = B[3]

(C[5, 2] + C[5, 6] + C[5, 7] + C[5, 8]) * T[5] 
                              - C[2, 5] * T[2] 
                              - C[7, 5] * T[7] 
                              - C[8, 5] * T[8] 
                              = C[7, 5] * T[7] + q[5] + h[5] * A[5] * (T[5] - Tref[5])
A[3, 0]*T[2] +                             + A[3, 4]*T[5] + A[3, 6]*T[7] + A[3, 7]*T[8] = B[4]

Next step: generate coefficient matrix as per equation [1] above. Compute the number of nodes with unknown temperature from node connectivity matrix. The coefficient matrix A[i, j] would be of size N x N where N is number of rows in node connectivity array defined above. phi_T = [row[0] for row in node_con_array] to get first column which contains node numbers for which temperature needs to be calculated.

Note:

• Arrays start with index zero but the lowest item in node connectivity matrix may be any number ≥ 1
• Coefficient matrix cannot have any element on diagonal row = zero
• The thermal resistance values are linked to node numbers and stored as rows and columns defined by connected nodes
• Thermal resistance matrix contains 0 values and hence a check is required to calculate the inverse of thermal resistances

  Get first column of each row as center node: row_i = node_con_mat[i]; n_num = row_i[0]

  Inner Loop Loops over second item onwards on each row

   Calculate diagonal item of coefficient matrix A[i, i]

 10.8 |  -4.0 |  -4.0 |  -0.8 |   0.0 |   0.0 |   0.0 | 
 -4.0 |   4.0 |   0.0 |   0.0 |   0.0 |   0.0 |   0.0 | 
 -4.0 |   0.0 |   4.0 |   0.0 |   0.0 |   0.0 |   0.0 | 
 -0.8 |   0.0 |   0.0 |  12.8 |  -4.0 |  -4.0 |   0.0 | 
  0.0 |   0.0 |   0.0 |  -4.0 | 3.995 |   0.0 |   0.0 | 
  0.0 |   0.0 |   0.0 |  -4.0 |   0.0 | 7.995 |  -4.0 | 
  0.0 |   0.0 |   0.0 |   0.0 |   0.0 |  -4.0 |   8.0 | 

Computed value of {b} excluding for those with known T:

Computed temperatures in [K]: for nodes in first column of node connectivity matrix above. Validation of result with hand calculation is in progress though values are very close to results produced by TNSolver (without HTC boundary conditions).

288.0
356.2
368.7
368.7
339.1
323.0
339.9
350.9
362.0
373.0

1. Give option to provide default length, cross-section area, material (thermal conductivity) for all resistors
2. Let user chose positive integers for numbering nodes and not enforcing start from 1 without any gap
3. Check correctness of inputs such as missing values, string instead of numbers...
4. Add radiative heat transfer from resistors
5. Add radiative heat transfer from nodes
6. Add convection from resistors
7. Add transient calculations
8. Allow comments
9. Single input file

The code can also be used to solve electrical networks where thermal conductivity 'k' needs to be replaced by electrical conductivity 'σ'. RTH = L/k/A, REL = ρ.L/A = L/σ/A where electrical conductivity σ = 1/ρ, ρ = electrical resistivity.

An iterative method like Jacobi or Gauss-Seidel iteration is required. The coefficient matrix A should ideally be diagonally dominant for guaranteed convergence. In Jacobi iteration, the value of each \(x_{i}\) at the new iteration is computed using the formula:

def Jacobi(A, b, x_0, max_iter=100, tolerance=1.0e-6):
  '''
  x_0: Initial guess for the solution, numpy array
  tolerance: Convergence criteria
  max_iter:  Maximum number of iterations to perform

  Returns: the approximate solution [vector].
  '''
  # Initialize the solution vector
  x = x_0.copy()
  n = len(b)
  for k in range(max_iter):
    x_new = x.copy()
    for i in range(n):
      # Calculate the sum of row elements, excluding diagonal value
      sum_row = 0
      for j in range(n):
        if i != j:
          sum_row =  sum_row + A[i, j] * x[j]
      x_new[i] = (b[i] - sum_row) / A[i, i]
    
    # Check for convergence, np.inf = Maximum absolute value
    L2NORM_new = np.linalg.norm(x_new - x, ord=np.inf)
    L2NORM_old = np.linalg.norm(x_new, ord=np.inf)
    if (L2NORM_new / L2NORM_old) <= tolerance:
      print(f"Solution converged after {k+1} iterations.")
      return x_new
      
    x = x_new
    
  print(f"Converge could not be achieved within {max_iter} iterations.")
  return x

Gauss-Seidel Method: The convergence is guaranteed if the matrix A is strictly diagonally dominant [magnitude of diagonal entry ≥ sum of the magnitudes of non-diagonal entries for each row] or symmetric positive definite [square matrix that is equal to its transpose]. The Gauss-Seidel method often converges faster than the Jacobi method by using updated values within the same iteration.

def GaussSeidel(A, b, x_0, max_iter=100, tolerance=1.0e-6, ):
  A = np.array(A, dtype=float)
  b = np.array(b, dtype=float)
  x = np.array(x_0, dtype=float)
  n = len(b)
  for iter in range(max_iter):
    x_old = np.copy(x)
    for i in range(n):
      sum_ax = np.dot(A[i, :i], x[:i]) + np.dot(A[i, i+1:], x_old[i+1:])
      x[i] = (b[i] - sum_ax) / A[i, i]
    if np.linalg.norm(x - x_old) < tolerance:
      print(f"Convergence achieved after {iter + 1} iterations.")
      return x
  print(f"Run did not converge within {max_iter} iterations.")
  return x

Finally, the most generic equation for thermal network solvers is derived below. Here, index 'i' is the node under iteration and 'j' refers nodes connected to it as described in node connectivity array.

\(\sum\limits_{j}\left[C_{ij} * (T_i - T_j)\right] + F_{ij} \cdot ε \cdot σ \cdot \sum\limits_{j}\left[\frac{A_{ij}}{2}) \cdot (T_{i}^4 - T_{REF}^4) \right] + \sum\limits_{j}\left[h \cdot \frac{A_{ij}}{2} * (T_i - T_j)\right] = \dot{q}_i \)

\(\sum\limits_{j}\left[ C_{ij} * (T_i - T_j)\right] + F_{ij} \cdot ε \cdot σ \cdot \sum\limits_{j}\left[\frac{A_{ij}}{2} \cdot (T_{i}^2 - T_{REF}^2) \cdot ((T_{i}^2 + T_{REF}^2)\right] + \sum\limits_{j}\left[h \cdot \frac{A_{ij}}{2} * (T_i - T_j)\right] = \dot{q}_i \)

\(\sum\limits_{j}\left[ C_{ij} * (T_i - T_j)\right] + F_{ij} \cdot ε \cdot σ \cdot \sum\limits_{j}\left[\frac{A_{ij}}{2} \cdot (T_{i} - T_{REF}) \cdot (T_{i} + T_{REF}) \cdot (T_{i}^2 + T_{REF}^2) \right] + \sum\limits_{j}\left[h \cdot \frac{A_{ij}}{2} * (T_i - T_j)\right] = \dot{q}_i\)

\(\sum\limits_{j}\left[C_{ij} + F_{ij} ε \cdotσ \cdot\frac{A_{ij}}{2}\cdot C_{RAD} + h \cdot \frac{A_{ij}}{2} \right ]\cdot T_i - \sum\limits_{j}\left (C_{ij} + h \cdot \frac{A{ij}}{2}) \right) \cdot T_j = \dot{q}_i + \sum\limits_{j}F_{ij} \cdot ε \cdot σ\cdot \frac{A_{ij}}{2} \left[T_{REF}\cdot C_{RAD} \right]\)

\(K_{ii}\cdot T_i - \sum\limits_{j}K_{ij}\cdot T_j = b_i \)

where \[K_{ii} = \sum\limits_{j}\left[C_{ij} + F_{ij} ε \cdotσ \cdot\frac{A_{ij}}{2}\cdot C_{RAD} + h \cdot \frac{A_{ij}}{2} \right ]\] \[ K_{ij} = - \left (C_{ij} + h \cdot\frac{A{ij}}{2} \right) \]

Sample thermal network for a battery pack

                   [CON1+]        [CON1-]          [CON2+]         [CON2-]
                     |              |                 |               |
                     |__[Case-top]__|                 |___[Case-top]__|    ...other cells     
                            |                                 |
                            |                                 |
                        [Cell-Z]                          [Cell-Z]
                            |                                 |
                            |                                 |
[Case]--[Gap-1]--[Can-1]--[C-1]--[Can-1]--[Gap-2]--[Can-2]--[C-2]--[Gap-3] ...other cells
                            |                                 |
                            |                                 |
                        [Cell-Z]                          [Cell-Z]
                            |                                 |
                            |                                 |
                        [Case-Bot]                        [Case-Bot]       ...other cells

TNSolver Examples

  |```````|
 [in]   [out]---[Tinf] h = 2.3 [W/m2-K]
  |       |
   ```````
Begin Solution Parameters
  title   = Simple Wall Model
  type    = steady
  units   = SI
  T units = C            !Default is 'C', other options K, F, R
End Solution Parameters

Begin Conductors
! label  type       nd_i   nd_j   parameters
  wall   conduction in     out    2.3  1.2  1.0   ! k L A
  fluid  convection out    Tinf   2.3  1.0        ! h A
End Conductors

Begin Boundary Conditions
! type      parameter   node(s)
  fixed_T   21.0        in                       ! Inner wall T
  fixed_T    5.0        Tinf                     ! Fluid T
End Boundary Conditions

 time = 0
    in  21
   out  12.2727
  Tinf  5

"label", "type", "nd_i", "nd_j", "T_i (C)", "T_j (C)", "Q (W)", "U (W/m^2-K)", "A (m^2)"
"wall",  "conduction", "in",  "out",  21, 12.2727, 16.7273, 1.91667, 1
"fluid", "convection", "out", "Tinf", 12.2727, 5,  16.7273, 2.3, 1

"label", "material", "volume (m^3)", "temperature (C)"
"in",    "N/A", 0, 21
"out",   "N/A", 0, 12.2727
"Tinf",  "N/A", 0, 5

Solution of the sample case described above (and used for Python scripting), in TNSolver without convection at nodes: note that TNSolver has option of specified temperature, heat flux, heat source and volumetric heat source at nodes. Convection can be specified on resistors only. Thus, tip convection can be applied using additional node to represent T_∞ or the reference temperature for HTC.

Begin Solution Parameters
  title   = Simple Wall Model
  type    = steady
  units   = SI
  T units = K            !Default is 'C', other options K, F, R
End Solution Parameters

Begin Conductors
! label  type        nd_i   nd_j   parameters
  R01    conduction  N_01   N_02   200  0.050  0.0005   ! k L A
  R02    conduction  N_02   N_03   200  0.025  0.0005   ! k L A
  R03    conduction  N_02   N_04   200  0.025  0.0005   ! k L A
  R04    conduction  N_02   N_05   200  0.125  0.0005   ! k L A
  R05    conduction  N_05   N_06   200  0.025  0.0005   ! k L A
  R06    conduction  N_05   N_07   200  0.025  0.0005   ! k L A
  R07    conduction  N_05   N_08   200  0.025  0.0005   ! k L A
  R08    conduction  N_08   N_09   200  0.025  0.0005   ! k L A
  R09    conduction  N_09   N_10   200  0.025  0.0005   ! k L A
End Conductors

Begin Boundary Conditions
! type      parameter   node(s)
  fixed_T   288         N_01
  fixed_T   323         N_06
  fixed_T   373         N_10
End Boundary Conditions

Begin Sources
  ! type  parameter(s)  node(s)
    Qsrc  50.0          N_02
    Qsrc  50.0          N_03
    Qsrc  50.0          N_04
End Sources

 time = 0
 N_01  288.0
 N_02  355.9
 N_03  368.4
 N_04  368.4
 N_05  338.2
 N_06  323.0
 N_07  338.2
 N_08  349.8
 N_09  361.4
 N_10  373.0

[288.0  355.9  368.4  368.4  338.2  323.0  338.2  349.8  361.4  373.0]

Example taken from "TNSolver: An Open Source Thermal Network Solver for Octave or MATLAB" by Bob Cochran at Thermal & Fluids Analysis Workshop (TFAWS) 2016

Begin Solution Parameters
  title = Radiation Heat Transfer Experiment - Black Target
  type  = steady
  nonlinear convergence = 1.0e-8
  maximum nonlinear iterations = 50
End Solution Parameters

Begin Conductors
  ! Conduction through the beaker wall, 0.1" thick pyrex glass
    t-bbin conduction targ bbin 0.14 0.00254 0.024829 ! k L A
  
  ! Convection from beaker to water
  ! label type        nd_i nd_j    mat     L       A
    t-w   ENChplateup bbin water   water   0.04445 0.02483
  
  ! Convection from target to air
  ! label type          nd_i   nd_j   mat  L       A
    t-air ENChplatedown targ   env    air  0.0889  0.0248

  ! Convection from outer shield to air
  ! label type          nd_i   nd_j   mat  H       L       theta  A
    s-air ENCiplateup   s_out  env    air  0.0508  0.0508  48.0   0.056439

  ! Conduction from inner to outer side of shield
    shield conduction s_in  s_out  steel  0.001   0.056439
End Conductors
  
Begin Radiation Enclosure
  ! surf   emiss  A        Fij
    htr    0.92   0.06701  0.0      0.1264   0.68415  0.0     0.1893
    targ   0.95   0.02482  0.34132  0.0      0.00603  0.1031  0.5494
    s_in   0.28   0.05643  0.81231  0.00265  0.12278  0.0     0.0622
    s_out  0.28   0.05643  0.0      0.0453   0.0      0.0     0.9546
    env    1.00   0.11840  0.10711  0.1151   0.02965  0.4547  0.2933
End Radiation Enclosure

Begin Boundary Conditions
  ! type     Tb     Node(s)
    fixed_T  23.0   env
    fixed_T  88.3   water
    fixed_T  515.0  htr
End Boundary Conditions

time = 0
bbin   93.4466
env    23
htr    515
s_in   272.032
s_out  271.956
targ   175.252
water  88.3

Fin with force convection and thermal radiation

                               [Tc]
                                |
            .,,,,,,,,,,,,,,,,,,,*,,,,,,,,,,,,,,,,,,,.      [Tc]
            |       |       |       |       |       |        |
            |       |       |       |       |       |        |
[base]-----[1]-----[2]-----[3]-----[4]-----[5]-----[6]-----[tip]
            |       |       |       |       |       |        |
            |       |       |       |       |       |        |
             ```````````````````|```````````````````       [Tr]
                              [Tr] 

! Fin Experiment Model - Aluminum fin - Forced convection

Begin Solution Parameters
  title = Aluminum Fin - Forced Convection
  type                         = steady
  nonlinear convergence        = 1.0e-8
  maximum nonlinear iterations = 15
End Solution Parameters

Begin Conductors
! label type       nd_i  nd_j   k        L        A
  100  conduction  base   1   167.0    0.030    0.000127 
  101  conduction   1     2   167.0    0.025    0.000127 
  102  conduction   2     3   167.0    0.045    0.000127 
  103  conduction   3     4   167.0    0.070    0.000127 
  104  conduction   4     5   167.0    0.050    0.000127 
  105  conduction   5     6   167.0    0.063    0.000127 
  106  conduction   6    tip  167.0    0.002    0.000127 

! label type      nd_i   nd_j mat  Vel   D        A
  121  EFCcyl      1     Tc   air  3.7  0.0127  0.00170
  122  EFCcyl      2     Tc   air  3.7  0.0127  0.00140
  123  EFCcyl      3     Tc   air  3.7  0.0127  0.00229
  124  EFCcyl      4     Tc   air  3.7  0.0127  0.00239
  125  EFCcyl      5     Tc   air  3.7  0.0127  0.00225
  126  EFCcyl      6     Tc   air  3.7  0.0127  0.00134
  127  EFCcyl     tip    Tc   air  3.7  0.0127  0.00013

! label type      nd_i   nd_j  emissivity     A
  221  surfrad     1     Tr      0.09    0.001696
  222  surfrad     2     Tr      0.09    0.001396
  223  surfrad     3     Tr      0.09    0.002294
  224  surfrad     4     Tr      0.09    0.002394
  225  surfrad     5     Tr      0.09    0.002254
  226  surfrad     6     Tr      0.09    0.001337
  227  surfrad    tip    Tr      0.09    0.000127
End Conductors

Begin Boundary Conditions
!  type        Tb      Node(s)
  fixed_T     25.0      Tc
  fixed_T     25.0      Tr
  fixed_T     55.8      base
End Boundary Conditions

Example of Composite Walls

     .---------,----------------------,-----------,
     |         |                      |           |
     |         |          kB          |           |
     |         |                      |           |
   [TIN]  kA  [1]--------------------[2]     kD  [TOUT]
     |         |                      |           |
     |         |          kC          |           |
     |_________|______________________|___________|

Option-1 of thermal network

             |````````|
             |        |
[TIN]-------[1]      [2]------[TOUT]
             |        |
             |        |
              ````````
! label type         node-1 node-2  parameters (k, L, A)
  100   conduction    Tin    1       1.0  1.0  2.0
  101   conduction    1      2       2.0  3.0  1.0
  102   conduction    1      2       2.0  3.0  1.0
  103   conduction    2      Tout    1.0  1.0  2.0


Option-2: Two parallel paths
  |`````````[1]``````[2]````````|
  |                             |
[TIN]                         [TOUT]
  |                             |
  |                             |
   `````````[3]``````[4]````````
! label type         node-1 node-2  parameters (k, L, A)
  100   conduction    Tin    1       1.0  1.0  1.0
  101   conduction    1      2       2.0  3.0  1.0
  102   conduction    2      Tout    1.0  1.0  1.0
  103   conduction    Tin    3       1.0  1.0  1.0
  104   conduction    3      4       2.0  3.0  1.0
  105   conduction    4      Tout    1.0  1.0  1.0

Fluid Network Solver

Major difference is in the basic definitions though: ΔT = Q. RTH and ΔP = Q² * RHD

1. Make an initial guess of flow rates in each flow resistors
2. Compute the sum of head losses around a loop of the network taking into account direction of flow: \(\sum \limits_{j}^{LOOP}\xi_j \cdot Q_{j}^n \)
3. Calculate the sum of absolute values of differential of head loss equation: \(\sum \limits_{j}^{LOOP}\lvert n\cdot \xi_j \cdot Q_{j}^{n-1} \rvert \)
4. Compute change in flow ΔQ by dividing results from step-2 by that of step-3.
5. Repeat steps 2 to 4 for each loop in the flow network
6. Repeat steps 2 to 5 iteratively until ΔQ computed during an iteration reduce below certain pre-defined tolerance (convergence criteria).

                [4]```````[5]``````[6]
                 |                  |
[1]-----[2]-----[3]                [7]
                 |                  |
                [8]-------[9]------[10]
                           |
                [12]------[11]-----[13]
                 |                  |
[17]-----------[14]                 |
                 |                  |
               [15]----------------[16]        

[FLOW RESISTORS BEGIN]
R_ID, s_node, e_node,  delta_p,  flow_rate, loss_coeff, Aref
R01,   1,       2,       0,         0,         0.02,    0.001
R02,   2,       3,       0,         0,         0.02,    0.001
R03,   3,       4,       0,         0,         0.02,    0.001
R04,   4,       5,       0,         0,         0.02,    0.001
R05,   5,       6,       0,         0,         0.02,    0.001
R06,   6,       7,       0,         0,         0.02,    0.001
R07,   7,      10,       0,         0,         0.02,    0.001
R08,   9,      10,       0,         0,         0.02,    0.001
R09,   8,       9,       0,         0,         0.02,    0.001
R10,   8,       3,       0,         0,         0.02,    0.001
R11,   9,      11,      20,         0,         0.02,    0.001
R12,  11,      12,       0,         0,         0.02,    0.001
R13,  11,      13,       0,      0.01,         0.02,    0.001
R14,  12,      14,       0,         0,         0.02,    0.001
R15,  14,      15,       0,         0,         0.02,    0.001
R16,  13,      16,       0,         0,         0.02,    0.001
R17,  15,      16,       0,         0,         0.02,    0.001
R18,  14,      17,       0,         0,         0.02,    0.001
[FLOW RESISTORS END]

[FLOW NODES BEGIN]
NNUM,  pressure,   pressure_loss
   1,       250,    0
  17,         0,    0
[FLOW NODES END]

             [m]         [n]
              |           |
              |           |
      Qi---> [i]---------[j] <---Qj
             /            \
            /              \
           [u]             [v]

Instead of tracking sign of flow rates (direction), it is easier to use following formulation:

Another approach is to use SIMPLE-like algorithm used in 3D CFD simulation programs. The "true" pressure and flow rates are defined as \(P = P^{*} + P^{\prime }\) and \(Q = Q^{*} + Q^{\prime }\). By linearizing the momentum equation, the flow correction \(Q_{ij}^{\prime }\) is related to the pressure correction \(P^{\prime }\) by: \(Q_{ij}^{\prime }\approx d_{ij}(P_{i}^{\prime } - P_{j}^{\prime })\) where \(d_{ij}\) is a sensitivity coefficient derived from the derivative of the head loss equation. The quadratic head loss (\(\Delta P = R\cdot Q^{2}\)) is linearized using the derivative \(d = \frac{1}{2R\cdot Q}\)

Step-1: Guess initial values of pressures at nodes

Step-3: Calculate flow rate correction Q' by differentiating branch pressure drop equation with respect to Q* and thus \(P'_i - P'_j = \frac {Q_{ij}}{|Q_{ij}|} \cdot \xi_{ij}\rho \cdot \frac {Q_{ij}^*}{A^2}\)

Step-4: Apply continuity equation at each node, excluding those at boundaries and known pressures to generate system of equation of the form [A]{P'} = {b}. Substitute \(Q_{ij} = Q_{ij}^{*} + Q_{ij}^{\prime }\) into the continuity equation at each node to obtain a system of linear equations for \(P^{\prime }\): \(\sum \limits_{j}d_{ij}(P_{i}^{\prime } - P_{j}^{\prime }) = S_{i} - \sum \limits_{j}Q_{ij}^{*}\). The right-hand side represents the mass imbalance (residual) at node 'i' and 1/β ≡ d_ij

Step-5: Solve for nodal values of P' and update nodal pressures, P = P* + P' and then compute element flow rate corrections and update branch flow rates: Q_ij = Q*_ij + Q'_ij

Step-6: Check for convergence of P and Q

Step-8: Go to step-2 if convergence not achieved.

Newton-Raphson Method for One Non-Linear Equation

For more explanations, refer to page engcourses-uofa.ca/books/numericalanalysis/nonlinear-systems-of-equations/newton-raphson-method

Example Validation Cases

Validation Case-1

[0]-------[1]-------[2]-------[3]
R_ID s_node e_node RHYD
R1   0      1      3
R2   1      2      5
R3   2      3      8

Node   Pressure
0      100
3      0

Validation Case-2

[0]-------[1]-------[2]-------[3]
                     |
                     |
                    [4]
R_ID s_node e_node RHYD
R1   0      1      3
R2   1      2      5
R3   2      3      8
R3   2      4      2

Node   Pressure
0      100
3      0
4      0

Step-1:

Step-2:

[[+1   0   0]
 [-1  +1   0]
 [ 0  -1  +1]
 [ 0   0  -1]]

Step-3:

Step-4:

0.5 / (R[0][1] * Q[0][1]) * (P0' - P1')             =  Q*01
0.5 / (R[1][2] * Q[1][2]) *       (P1' - P2')       =  Q*12
0.5 / (R[2][3] * Q[2][3]) *             (P2' - P3') =  Q*23
or
1/β01|+1  -1   0   0| |P0'|    |b0|
1/β12| 0  +1  -1   0| |P1'| =  |b1|
1/β23| 0   0  +1  -1| |P2'|    |b2|
                      |P3'|

[A] =

[[  1       0          0        0   ]
 [  0       100.03   -100.0     0   ]
 [  0      -100.00    162.5     0   ]
 [  0.      0          0        1   ]]

Step-5:

Both the node pressures and flow rates though branches match the hand calculations for the "Validation Case-1". You may test the following code for your application and let me know cases where it fails.

import numpy as np
def tupleToList(list_of_tuples):
  '''
  Return the length of unique items in a tuple defining branches
  '''
  list_1, list_2 = [], []
  for t in list_of_tuples:
    list_1.append(t[0])
    list_2.append(t[1])
  unique_list = list(set(list(set(list_1)) + list(set(list_2))))
  unique_list.sort()
  return len(unique_list), list_1, list_2
  
def pipeNetwork(pipes, pr_bc, alpha_p=0.5, alpha_q=0.7, max_iter=500, tol=1e-5):
  """
  SIMPLE algorithm for arbitrary pipe networks with non-linear resistances 
  and multiple pressure boundary conditions. Parameters are:
    pipes: list of (m, n, R) - Pipe from m -> n with resistance R
    pr_bc: dictionary - {node_index: pressure_value}
  """
  eps = 1e-8
  Np = len(pipes)
  nn_nodes, ni, nj = tupleToList(pipes)
  nn_uknown = nn_nodes - len(pr_bc.items())

  # Incidence matrix
  B = np.zeros((nn_nodes, Np))  # N_nodes x N_branch = rows x cols
  R = np.zeros(Np)

  # Initial guess for pressures at nodes
  p = np.ones(nn_nodes)
  for node, val in pr_bc.items():
    p[node] = val
  Q = np.ones(Np) * 0.001   # small non-zero initial flow
  
  for i, (m, n, Ri) in enumerate(pipes):
    B[m, i] = +1.0
    B[n, i] = -1.0
    R[i] = Ri
    if p[m] - p[n] != 0:
      Q[i] = np.sign(p[m] - p[n]) * np.sqrt(abs(p[m] - p[n]) / R[i])

  # Start SIMPLE-type iterations
  for it in range(max_iter):

    # 1. Update effective resistance (non-linearity)
    R_eff = R * np.abs(Q) + eps
    
    # 2. Momentum predictor
    dp = B.T @ p
  
    Q_star = dp / R_eff  # Q * Q = dp/R, Q = dp / [R * Q] = dp/R_eff

    # 3. Mass imbalance
    b = B @ Q_star

    # 4. Pressure correction system: dp/dq = 2RQ = 0.5/R_eff
    d = 0.5 / R_eff
    Ap = B @ np.diag(d) @ B.T
    
    # Apply pressure boundary conditions
    for node, val in pr_bc.items():
      Ap[node, :] = 0.0
      Ap[:, node] = 0.0
      Ap[node, node] = 1.0
      b[node] = 0.0
   
    # 5. Solve pressure correction: "source term" is the mass imbalance
    p_corr = np.linalg.solve(Ap, -b)
    
    # 6. Update pressures
    p = p + alpha_p * p_corr
    for node, val in pr_bc.items():
      p[node] = val
    
    # 7. Correct flows
    for i, (m, n, Ri) in enumerate(pipes):
      if p[m] - p[n] != 0:
        Q[i] = np.sign(p[m] - p[n]) * np.sqrt(abs(p[m] - p[n]) / R[i])
    
    # 8. Check convergence: return node connectivity, node pressure, flow rates
    if np.max(np.abs(b)) < tol:
      print(f"Converged in {it+1} iterations")
      return ni, nj, p, Q
  
  # Print divergence and return 'None'
  print("SIMPLE did not converge!")
  return ni, nj, None, None
  
pipes = [(0, 1, 3.0), (1, 2, 5.0), (2, 3, 8.0), (2, 4, 2.0)]
pr_bc = {0: 100.0, 3: 0.0, 4: 0.0}

ni, nj, nodeP, Qdot = pipeNetwork(pipes=pipes, pr_bc=pr_bc, max_iter=100)
if nodeP is not None:
  print("-" * 50)
  print("All values are in chosen unit of pressure and [m^3/s]")
  print("Node pressures:")
  for index, value in enumerate(nodeP, start=1):
    print(f"{index}: {value:8.2f}")

if Qdot is not None:
  print("-" * 50)
  print("Flow through pipe branches:")
  for index, value in enumerate(Qdot):
    print(f" {ni[index]} -> {nj[index]} = {value:6.3f}  |", end="")
print("\n")

Hardy-Cross Method Example:

Step 1: Assign Initial Flow Rates (Q) and R values

Pipe   Length (L)   Diameter (D)   Friction Factor (f)   Resistance (R)  Flow (Q)
  1    ...             ...              ...              0.5             10
  2    ...             ...              ...              0.8             6
  3    ...             ...              ...              1.2             4

Step 2: Calculate the Head Loss (\(\mathbf{h}_{\mathbf{f}}\)) for each pipe The head loss is \(h_{f}=R\cdot Q|Q|\) to maintain the correct sign for head loss regardless of flow direction. Pipe 1: \(h_{f1}=R_{1}\cdot Q_{1}|Q_{1}|=0.5\cdot 10\cdot 10=+50\), Pipe 2: \(h_{f2}=R_{2}\cdot Q_{2}|Q_{2}|=0.8\cdot 6\cdot 6=+28.8\), Pipe 3: \(h_{f3}=R_{3}\cdot Q_{3}|Q_{3}|=1.2\cdot (-4)\cdot |-4|=-19.2\) (Negative because flow is counter-clockwise)

Step 3: Check the Head Loss Summation around the loop According to energy conservation, the sum of head losses around a closed loop should be zero (\(\sum h_{f}=0\)). \(\sum h_{f}=h_{f1}+h_{f2}+h_{f3}=50+28.8+(-19.2)=59.6\) The result is not zero, so the initial flow rates are incorrect.

Step 5: Update the Flow Rates Apply the flow correction to each pipe in the loop. The sign of \(\Delta Q\) is negative, meaning the flows in the clockwise direction were overestimated and need to be reduced, or the counter-clockwise flow needs to be increased. Pipe 1 (Clockwise): \(Q_{1,new}=Q_{1}+\Delta Q=10+(-2.04)=7.96\), Pipe 2 (Clockwise): \(Q_{2,new} = Q_{2}+\Delta Q = 6+(-2.04)=3.96\), Pipe 3 (Counter-clockwise, so subtract \(\Delta Q\)): \(Q_{3,new}=Q_{3}-\Delta Q=-4-(-2.04)=-1.96\)

Example: 3-Reservoirs System, here \(R=\frac{8fL}{\pi^{2}gD^{5}}\)

|       |   
|  [1]  |             |       |
|_______|             |  [3]  | 
   \           /``````|_______|
    \         /
     \       /
      \     /
       \   /             |       |
        \ /              |  [2]  |
      [4]*---------------|_______|
          \
           \
===========---=====================Ground==========
H = (70   80  100   ?) [m], friction factor = 0.015
L[1, 4] =  5000 [m], L[2, 4] = 3000 [m], L[3, 4] = 4000 [m]
D[1, 4] =  0.60 [m], D[2, 4] = 0.80 [m], D[3, 4] = 1.20 [m]
R[1, 4] = 72.7266,   R[2, 4] = 11.3517 , R[3, 4] = 1.99316

Pipe [1, 4]: H[1] - H[4] = R[1, 4] * Q[1, 4] * |Q[1, 4]|
Pipe [2, 4]: H[2] - H[4] = R[2, 4] * Q[2, 4] * |Q[2, 4]|
Pipe [3, 4]: H[3] - H[4] = R[3, 4] * Q[3, 4] * |Q[3, 4]|

Pipe [1, 4]:  70 = H[4] + 72.7266 * Q[1, 4] * |Q[1, 4]|
Pipe [2, 4]:  80 = H[4] + 1.99316 * Q[2, 4] * |Q[2, 4]|
Pipe [3, 4]: 100 = H[4] + 11.3517 * Q[3, 4] * |Q[3, 4]|

H[4]    R[1, 4]*Q[1, 4]*|Q[1, 4]|               0                           0               = 70
H[4]                0               R[2, 4]*Q[2, 4]*|Q[2, 4]|                           0   = 100
H[4]                0                           0               R[3, 4]*Q[3, 4]*|Q[3, 4]|   = 80
   0            Q[1, 4]                     Q[2, 4]                     Q[3, 4]             = 0

[1    R[1, 4]*|Q[1, 4]|                    0                     ] {H[4]   }    {H[1]}   {e1}
[1                    0    R[2, 4]*|Q[2, 4]|                     ] {Q[1, 4]} _  {H[2]} = {e2}
[1                    0                    0    R[3, 4]*|Q[3, 4]|] {Q[2, 4]}    {H[3]}   {e3}
[0                    1                    1                    1] {Q[3, 4]}    {   0]   {e4}

Results are: Q[1, 4] = -0.395 [m³/s] i.e. 0.395 [m³/s] from node 4 to reservoir 1 (opposite the original assumed direction), Q[2, 4] = +1.273 [m³/s] from reservoir 2 to the node 4 (along the assumed flow direction in the pipe), Q[3, 4] = -0.890 [m³/s] i.e. 0.890 [m³/s] from the node 4 to reservoir 3 (opposite the original assumed direction) and h[4] = 81.60 [m] with mass imbalance of -0.012 [m³/s] at node '4'.

pipes = [(0, 1, 72.7266), (0, 2, 11.3517), (0, 3, 1.99316)]
pr_bc = {1: 70.0, 2: 100.0, 3: 80.0}

Converged in 18 iterations
--------------------------------------------------
All values are in chosen unit of pressure and [m^3/s]
Node pressures:
1:    81.53
2:    70.00
3:   100.00
4:    80.00
--------------------------------------------------
Flow through pipe branches:
 0 -> 1 =  0.398  | 0 -> 2 = -1.275  | 0 -> 3 =  0.877  |

References:

• A New Approach to Modeling Fluid/Gas Flows in Networks by W. S. Winters, SANDIA REPORT SAND2001-8422
• Algorithms for Pipe Network Analysis and Their Reliability by Don J. Wood, University of Kentucky
• PIPE NETWORK ANALYSIS By Mun-Fong Lee, University of Florida, Gainesville