add friction model "Hofer"

CHANGELOG.rst

@@ -1,6 +1,13 @@
 Change Log
 =============
 
+[upcoming release] - 2025-xx-xx
+-------------------------------
+- [ADDED] Hofer formula for friction factor "lambda"
+- [ADDED] deprecation warning for "max_iter" (use "max_iter_hyd" or "max_iter_therm" instead)
+- [CHANGED] default maximum iterations for hydraulic calculations now 30 instead of 10
+
+
 [0.11.0] - 2024-11-07
 -------------------------------
 - [ADDED] heat_consumer plotting
@@ -28,8 +35,6 @@ Change Log
 - [FIXED] alpha also applied to mdot
 - [REMOVED] support for Python 3.8 due to EOL
 
-
-
 [0.10.0] - 2024-04-09
 -------------------------------
 

doc/source/components/pipe/pipe_component.rst

@@ -142,46 +142,65 @@ pipe sections.
 Friction models
 ^^^^^^^^^^^^^^^
 
-Three friction models are used to calculate the velocity dependent friction factor:
+Four friction models are available to calculate the velocity dependent friction factor :math:`\lambda`:
 
-- Nikuradse
-- Prandtl-Colebrook
-- Swamee-Jain
+- Nikuradse ("nikuradse")
+- Prandtl-Colebrook ("colebrook")
+- Hofer ("hofer")
+- Swamee-Jain ("swamee-jain")
 
-Nikuradse is chosen by default. In this case, the friction factor is calculated by:
+They are set by the :code:`friction_model` parameter and the name given in parentheses.
+*Nikuradse* is chosen by default. In this case, the friction factor is calculated by:
 
 .. math::
    :nowrap:
 
    \begin{align*}
-    \lambda &= \frac{64}{Re} + \frac{1}{(-2 \cdot \log (\frac{k}{3.71 \cdot d}))^2}\\
+    \lambda &= \frac{64}{Re} + \frac{1}{\left(-2 \cdot \log \left(\frac{k}{3.71 \cdot d}\right)\right)^2}\\
    \end{align*}
 
 
 Note that in literature, Nikuradse is known as a model for turbulent flows. In pandapipes, the formula for the
 Nikuradse model is also applied for laminar flow.
 
-If Prandtl-Colebrook is selected, the friction factor is calculated iteratively according to
+If *Prandtl-Colebrook* (also known as Colebrook-White) is selected, the friction factor is calculated iteratively according to
 
 .. math::
    :nowrap:
 
    \begin{align*}
-     \frac{1}{\sqrt{\lambda}} &= -2 \cdot \log (\frac{2.51}{Re \cdot \sqrt{\lambda}} + \frac{k}{3.71 \cdot d})\\
+     \frac{1}{\sqrt{\lambda}} &= -2 \cdot \log \left(\frac{2.51}{Re \cdot \sqrt{\lambda}} + \frac{k}{3.71 \cdot d}\right)\\
    \end{align*}
 
 Equations for pressure losses due to friction were taken from :cite:`Eberhard1990` and
 :cite:`Cerbe2008`.
 
-The equation according to Swamee-Jain :cite:`Swamee1976` is an approximation of the calculation method according
+The *Hofer-Equation* is an explicit approximation of the Prandtl-Colebrook method and defined as shown in
+:cite:`Benner.2019` (based on :cite:`Hofer.1973`)
+
+.. math::
+   :nowrap:
+
+   \begin{align*}
+    \lambda_H &= \frac{1}{(-2 \log \left(\frac{4.518}{Re} \cdot \log \left(\frac{Re}{7}) + \frac{k}{3.71 \cdot d}\right)\right)^2}\\
+   \end{align*}
+
+Very small Reynolds numbers can lead to problems with the two nested log-functions.
+Thus, in pandapipes, the laminar :math:`\lambda_l = 64/Re` is applied for small Reynolds numbers
+(:math:`Re < \underline{Re}`) and :math:`\lambda_H` is used for high Reynolds numbers (:math:`Re > \overline{Re}`).
+In the transition area :math:`\underline{Re} \leq Re \leq \overline{Re}`, :math:`\lambda_H` is interpolated linearly.
+By default, the boundaries of this transition area are set to :math:`\underline{Re}=2000` and :math:`\overline{Re}=3000`.
+Note that the derivative used for Hofer is the same as for Prandtl-Colebrook.
+
+The equation according to *Swamee-Jain* :cite:`Swamee1976` is another approximation of the calculation method according
 to Prandtl-Colebrook. It is an explicit formula for the friction factor of the transition
 zone of turbulent flows in pipes and is defined as follows:
 
 .. math::
    :nowrap:
 
    \begin{align*}
-    \lambda &= \frac{0.25}{(\log(\frac{k}{3.7 \cdot d} + \frac{5.74}{Re^{0.9}}))^2}\\
+    \lambda &= \frac{0.25}{\left(\log\left(\frac{k}{3.7 \cdot d} + \frac{5.74}{Re^{0.9}}\right)\right)^2}\\
    \end{align*}
 
 

doc/source/pipeflow/calculation_modes.rst

@@ -92,8 +92,8 @@ In gas flows, the velocity is typically not constant along a pipeline. For this
 tables for pipes show more entries in comparison with the result tables for incompressible media.
 
 
-Temperature calculations (pipeflow option: mode = "sequential", mode = "bidrectional" or mode = "heat")
-=======================================================================================================
+Temperature calculations (pipeflow option: mode = "sequential", mode = "bidirectional" or mode = "heat")
+========================================================================================================
 
 Important parameters of the network main components (junctions and pipes) needed for the calculation
 are listed in the following table. The :ref:`component section <components>` of this manual contains

doc/source/references.bib

@@ -81,4 +81,40 @@ @article{Schmidt.2015
  issn = {1389-4420},
  journal = {Optimization and Engineering},
  doi = {10.1007/s11081-014-9246-x}
-}
+}
+
+
+@incollection{Benner.2019,
+ author = {Benner, Peter and Grundel, Sara and Himpe, Christian and Huck, Christoph and Streubel, Tom and Tischendorf, Caren},
+ title = {Gas Network Benchmark Models},
+ pages = {171--197},
+ publisher = {{Springer International Publishing}},
+ isbn = {978-3-030-03717-8},
+ series = {Differential-Algebraic Equations Forum},
+ editor = {Campbell, Stephen and Ilchmann, Achim and Mehrmann, Volker and Reis, Timo},
+ booktitle = {Applications of Differential-Algebraic Equations: Examples and Benchmarks},
+ year = {2019},
+ address = {Cham},
+ doi = {10.1007/11221_2018_5},
+}
+
+@book{Campbell.2019,
+ year = {2019},
+ title = {Applications of Differential-Algebraic Equations: Examples and Benchmarks},
+ address = {Cham},
+ publisher = {{Springer International Publishing}},
+ isbn = {978-3-030-03717-8},
+ series = {Differential-Algebraic Equations Forum},
+ editor = {Campbell, Stephen and Ilchmann, Achim and Mehrmann, Volker and Reis, Timo},
+ doi = {10.1007/978-3-030-03718-5}
+}
+
+@article{Hofer.1973,
+ author = {Hofer, P.},
+ year = {1973},
+ title = {{Beurteilung von Fehlern in Rohrnetzberechnungen}},
+ pages = {113--119},
+ volume = {114},
+ number = {3},
+ journal = {Das Gas- und Wasserfach. GWF - Ausgabe Wasser, Abwasser}
+}

src/pandapipes/pf/derivative_calculation.py

@@ -140,26 +140,40 @@ def get_derived_values(node_pit, from_nodes, to_nodes, use_numba):
     return calc_derived_values_np(node_pit, from_nodes, to_nodes)
 
 
+def calc_lambda_transition_area(lambda_laminar, lambda_turb, re, begin_transition_re=2000,
+                                end_transition_re=3000):
+    """
+    Linear interpolation in the transition area between laminar and turbulent flow. Critical
+    point is usually at Reynolds = 2300, thus the default transition area is around that point.
+    """
+    lambda_tot = lambda_laminar
+    share_turb = (re - begin_transition_re) / (end_transition_re - begin_transition_re)
+    lambda_transition = (1 - share_turb) * lambda_laminar + share_turb * lambda_turb
+    lambda_tot[(re >= begin_transition_re) & (re < end_transition_re)] = \
+        lambda_transition[(re >= begin_transition_re) & (re < end_transition_re)]
+    lambda_tot[re >= end_transition_re] = lambda_turb[re >= end_transition_re]
+    return lambda_tot
+
 def calc_lambda(m, eta, d, k, gas_mode, friction_model, lengths, options, area):
     """
-    Function calculates the friction factor of a pipe. Turbulence is calculated based on
-    Nikuradse. If v equals 0, a value of 0.001 is used in order to avoid division by zero.
-    This should not be a problem as the pressure loss term will equal zero (lambda * u^2).
+    Function calculates the friction factor of a pipe, considering laminar and turbulent flow.
+    Turbulence is calculated based on the friction model provided. If m equals 0, a value of 0.001 is used in order
+    to avoid division by zero. This should not be a problem as the pressure loss term will equal zero (lambda * u^2).
 
-    :param m:
+    :param m: mass flow rate
     :type m:
     :param eta:
     :type eta:
     :param rho:
     :type rho:
-    :param d:
+    :param d: inner pipe diameter
     :type d:
-    :param k:
+    :param k: roughness coefficient
     :type k:
     :param gas_mode:
     :type gas_mode:
-    :param friction_model:
-    :type friction_model:
+    :param friction_model: friction formula, one of "nikuradse", "colebrook", "hofer", "swamee-jain"
+    :type friction_model: str
     :param lengths:
     :type lengths:
     :param options:
@@ -170,17 +184,19 @@ def calc_lambda(m, eta, d, k, gas_mode, friction_model, lengths, options, area):
     if options["use_numba"]:
         from pandapipes.pf.derivative_toolbox_numba import calc_lambda_nikuradse_incomp_numba as \
             calc_lambda_nikuradse_incomp, colebrook_numba as colebrook, \
-            calc_lambda_nikuradse_comp_numba as calc_lambda_nikuradse_comp
+            calc_lambda_nikuradse_comp_numba as calc_lambda_nikuradse_comp, \
+            calc_lambda_hofer_comp_numba as calc_lambda_hofer_comp
     else:
         from pandapipes.pf.derivative_toolbox import calc_lambda_nikuradse_incomp_np as \
             calc_lambda_nikuradse_incomp, colebrook_np as colebrook, \
-            calc_lambda_nikuradse_comp_np as calc_lambda_nikuradse_comp
+            calc_lambda_nikuradse_comp_np as calc_lambda_nikuradse_comp, \
+            calc_lambda_hofer_comp_np as calc_lambda_hofer_comp
     if gas_mode:
         re, lambda_laminar, lambda_nikuradse = calc_lambda_nikuradse_comp(m, d, k, eta, area)
     else:
         re, lambda_laminar, lambda_nikuradse = calc_lambda_nikuradse_incomp(m, d, k, eta, area)
 
-    if friction_model == "colebrook":
+    if friction_model.lower() == "colebrook":
         # TODO: move this import to top level if possible
         from pandapipes.pipeflow import PipeflowNotConverged
         max_iter = options.get("max_iter_colebrook", 100)
@@ -192,9 +208,14 @@ def calc_lambda(m, eta, d, k, gas_mode, friction_model, lengths, options, area):
                 "inconsistencies. The maximum iterations can be given as 'max_iter_colebrook' "
                 "argument to the pipeflow.")
         return lambda_colebrook, re
-    elif friction_model == "swamee-jain":
+    elif friction_model.lower() == "swamee-jain":
         lambda_swamee_jain = 0.25 / ((np.log10(k / (3.7 * d) + 5.74 / (re ** 0.9))) ** 2)
         return lambda_swamee_jain, re
+    elif friction_model.lower() == "hofer":
+        re, lambda_laminar, lambda_hofer = calc_lambda_hofer_comp(m, d, k, eta, area)
+        lambda_tot = calc_lambda_transition_area(lambda_laminar, lambda_hofer, re,
+                                                 begin_transition_re=2000, end_transition_re=4000)
+        return lambda_tot, re
     else:
         # lambda_tot = np.where(re > 2300, lambda_laminar + lambda_nikuradse, lambda_laminar)
         lambda_tot = lambda_laminar + lambda_nikuradse
@@ -203,12 +224,11 @@ def calc_lambda(m, eta, d, k, gas_mode, friction_model, lengths, options, area):
 
 def calc_der_lambda(m, eta, d, k, friction_model, lambda_pipe, area):
     """
-    Function calculates the derivative of lambda with respect to v. Turbulence is calculated based
-    on Nikuradse. This should not be a problem as the pressure loss term will equal zero
-    (lambda * u^2).
+    Function calculates the derivative of lambda with respect to m (mass flow rate). Turbulence is
+    calculated based on the friction model provided.
 
-    :param v:
-    :type v:
+    :param m:
+    :type m:
     :param eta:
     :type eta:
     :param rho:
@@ -229,9 +249,10 @@ def calc_der_lambda(m, eta, d, k, friction_model, lambda_pipe, area):
     df_dm = np.zeros_like(m)
     df_dlambda = np.zeros_like(m)
     lambda_der = np.zeros_like(m)
-    pos = m != 0
+    pos = ~np.isclose(abs(m), 0)
 
-    if friction_model == "colebrook":
+    if friction_model.lower() in ["colebrook", "hofer"]: # hofer is explicit approximation of
+        # colebrook
         b_term[pos] = (2.51 * eta[pos] * area[pos] / (m[pos] * d[pos] * np.sqrt(lambda_pipe[pos])) +
                        k[pos] / (3.71 * d[pos]))
 
@@ -244,7 +265,7 @@ def calc_der_lambda(m, eta, d, k, friction_model, lambda_pipe, area):
         lambda_der[pos] = df_dm[pos] / df_dlambda[pos]
 
         return lambda_der
-    elif friction_model == "swamee-jain":
+    elif friction_model.lower() == "swamee-jain":
         param = (k[pos] / (3.7 * d[pos]) + 5.74 * ((eta[pos] * area[pos]) /
                  (np.abs(m[pos]) * d[pos])) ** 0.9)
         # 0.5 / (log(10) * log(param)^3 * param) * 5.166 * abs(eta)^0.9  / (abs(rho * d)^0.9

src/pandapipes/pf/derivative_toolbox.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 from numpy import linalg
+from copy import deepcopy
 
 from pandapipes.constants import P_CONVERSION, GRAVITATION_CONSTANT, NORMAL_PRESSURE, \
     NORMAL_TEMPERATURE
@@ -73,24 +74,45 @@ def derivatives_hydraulic_comp_np(node_pit, branch_pit, lambda_, der_lambda, p_i
     return load_vec, load_vec_nodes_from, load_vec_nodes_to, df_dm, df_dm_nodes, df_dp, df_dp1
 
 
-def calc_lambda_nikuradse_incomp_np(m, d, k, eta, area):
+def calc_lambda_laminar(m, d, eta, area):
     m_abs = np.abs(m)
     re = np.divide(m_abs * d, eta * area)
     lambda_laminar = np.zeros_like(m)
     lambda_laminar[~np.isclose(re, 0)] = 64 / re[~np.isclose(re, 0)]
+    return lambda_laminar, re
+
+
+def calc_lambda_nikuradse_incomp_np(m, d, k, eta, area):
+    lambda_laminar, re = calc_lambda_laminar(m, d, eta, area)
     lambda_nikuradse = np.divide(1, (-2 * np.log10(k / (3.71 * d))) ** 2)
     return re, lambda_laminar, lambda_nikuradse
 
 
 def calc_lambda_nikuradse_comp_np(m, d, k, eta, area):
-    m_abs = np.abs(m)
-    re = np.divide(m_abs * d, eta * area)
-    lambda_laminar = np.zeros_like(m)
-    lambda_laminar[~np.isclose(re, 0)] = 64 / re[~np.isclose(re, 0)]
+    lambda_laminar, re = calc_lambda_laminar(m, d, eta, area)
     lambda_nikuradse = np.divide(1, (2 * np.log10(d / k) + 1.14) ** 2)
     return re, lambda_laminar, lambda_nikuradse
 
 
+def calc_lambda_hofer_comp_np(m, d, k, eta, area, hofer_re_threshold=2000):
+    """Calculate Lambda acc. to Hofer-Equation (explicit version of Colebrook) [1].
+    Due to nested log10 operations, small Reynolds numbers must be avoided.
+    It is justified because the laminar Lambda (instead of Hofer) will be applied for small
+    Reynolds numbers anyway.
+
+    [1]: P. Hofer. Beurteilung von Fehlern in Rohrnetzberechnungen (Error evaluation in calculation
+        of pipelines). GWF–Gas/Erdgas, 114(3):113–119, 1973
+    """
+    lambda_laminar, re = calc_lambda_laminar(m, d, eta, area)
+    lambda_hofer = np.zeros_like(re)
+    re_lower_lim = deepcopy(re)
+    re_lower_lim[re < hofer_re_threshold] = hofer_re_threshold
+    lambda_hofer = \
+        np.divide(1, (-2 * np.log10((4.518/re_lower_lim) * np.log10(re_lower_lim/7) + (k / (3.71 * d)))) ** 2)
+    lambda_hofer[re < hofer_re_threshold] = 0
+    return re, lambda_laminar, lambda_hofer
+
+
 def calc_medium_pressure_with_derivative_np(p_init_i_abs, p_init_i1_abs):
     val = 2 / 3
     p_m = p_init_i_abs.copy()

src/pandapipes/pf/derivative_toolbox_numba.py

@@ -111,15 +111,41 @@ def calc_lambda_nikuradse_comp_numba(m, d, k, eta, area):
     lambda_nikuradse = np.empty_like(m)
     lambda_laminar = np.zeros_like(m)
     re = np.empty_like(m)
+    m_abs = np.abs(m)
     for i, mi in enumerate(m):
-        m_abs = np.abs(mi)
-        re[i] = np.divide(m_abs * d[i], eta[i] * area[i])
+        re[i] = (m_abs[i] * d[i]) / (eta[i] * area[i])
         if re[i] != 0:
-            lambda_laminar[i] = np.divide(64, re[i])
-        lambda_nikuradse[i] = np.divide(1, (2 * np.log10(np.divide(d[i], k[i])) + 1.14) ** 2)
+            lambda_laminar[i] = 64 / re[i]
+        lambda_nikuradse[i] = (2 * np.log10((d[i] / k[i])) + 1.14) ** -2
     return re, lambda_laminar, lambda_nikuradse
 
 
+@jit((float64[:], float64[:], float64[:], float64[:], float64[:]), nopython=True)
+def calc_lambda_hofer_comp_numba(m, d, k, eta, area):
+    """Calculate Lambda acc. to Hofer-Equation (explicit version of Colebrook) [1].
+       Due to nested log10 operations, small Reynolds numbers must be avoided.
+       It is justified because the laminar Lambda (instead of Hofer) will be applied for small
+       Reynolds numbers anyway.
+
+       [1]: P. Hofer. Beurteilung von Fehlern in Rohrnetzberechnungen (Error evaluation in calculation
+           of pipelines). GWF–Gas/Erdgas, 114(3):113–119, 1973
+       """
+    hofer_re_threshold = 2000
+    lambda_hofer = np.zeros_like(m)
+    lambda_laminar = np.zeros_like(m)
+    re = np.empty_like(m)
+    m_abs = np.abs(m)
+    for i, mi in enumerate(m):
+        re[i] = (m_abs[i] * d[i]) / (eta[i] * area[i])
+        if re[i] != 0:
+            lambda_laminar[i] = 64 / re[i]
+        if re[i] >= hofer_re_threshold:
+            lambda_hofer[i] = (-2 * np.log10((4.518/re[i]) * np.log10(re[i]/7)
+                                             + (k[i] / (3.71 * d[i])))
+                               ) ** -2
+    return re, lambda_laminar, lambda_hofer
+
+
 @jit((float64[:], float64[:]), nopython=True, cache=False)
 def calc_medium_pressure_with_derivative_numba(p_init_i_abs, p_init_i1_abs):
     p_m = p_init_i_abs.copy()