Introduce FrictionFactorModel protocol with built‑in implementations and user extensibility

CHANGELOG.rst

@@ -2,6 +2,15 @@ Change Log
 =============
 [upcoming release] - 2026-..-..
 -------------------------------
+- [ADDED]
+  FrictionFactorModel protocol and RegimeAwareFrictionFactorModel.
+  Users can now implement their own friction factor models and pass them
+  to the pipeflow.  A regime‑aware model allows combining different
+  models for laminar, transient, and turbulent flow.
+
+- [CHANGED]
+  Friction model selection is now object‑based.  Passing a string
+  name (e.g. "colebrook") is still supported for backward compatibility.
 
 [0.14.0] - 2026-05-26
 -------------------------------

doc/source/components/pipe/pipe_component.rst

@@ -139,50 +139,13 @@ pipe sections.
 
 
 
-Friction models
-^^^^^^^^^^^^^^^
+Friction factor models
+^^^^^^^^^^^^^^^^^^^^^^
 
-Three friction models are used to calculate the velocity dependent friction factor:
-
-- Nikuradse
-- Prandtl-Colebrook
-- Swamee-Jain
-
-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}\\
-   \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
-
-.. math::
-   :nowrap:
-
-   \begin{align*}
-     \frac{1}{\sqrt{\lambda}} &= -2 \cdot \log (\frac{2.51}{Re \cdot \sqrt{\lambda}} + \frac{k}{3.71 \cdot d})\\
-   \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
-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}\\
-   \end{align*}
+The pipeflow solver uses a friction factor model to compute the Darcy‑Weisbach
+friction factor :math:`\lambda` and its derivative with respect to mass flow.
+Several built‑in models are available; you can also provide a custom model.
+For details, see :ref:`friction_factor_models`.
 
 
 Heat transfer mode

doc/source/conf.py

@@ -251,3 +251,10 @@
 
 # Additional magic lines (beacuse of: https://github.com/phn/pytpm/issues/3)
 numpydoc_show_class_members = False
+
+# needed to preserve typealiases in the docs
+autodoc_type_aliases = {
+    "Float64_1D": "Float64_1D",
+    "FrictionFactorResult": "FrictionFactorResult",
+    "LambdaEstimator": "LambdaEstimator",
+}

doc/source/pipeflow.rst

@@ -17,11 +17,10 @@ are used and how it is possible to influence the calculation.
 
 .. toctree::
     :maxdepth: 1
-    
+
     pipeflow/run
     pipeflow/options
     pipeflow/pipeflow_procedure
+    pipeflow/friction_factor_models
     pipeflow/calculation_modes
     pipeflow/internal_functions
-
-

doc/source/pipeflow/friction_factor_models.rst

@@ -0,0 +1,65 @@
+.. currentmodule:: pandapipes.pf.friction_factor_model
+
+.. _friction_factor_models:
+
+**********************
+Friction Factor Models
+**********************
+
+The friction factor :math:`\lambda` appears in the Darcy‑Weisbach type
+pressure drop equations used by pandapipes. The exact form of the pressure
+drop equation depends on the fluid model (incompressible or compressible).
+
+All friction factor models implement the :class:`~FrictionFactorModel`
+protocol and provide both :math:`\lambda` and its derivative
+:math:`\frac{\mathrm{d}\lambda}{\mathrm{d}\dot{m}}`,
+which is essential for :ref:`constructing the Jacobian matrix <jacobian>`
+in the Newton‑Raphson solver.
+
+.. note::
+    Because the friction factor depends only on the magnitude of the mass flow
+    (via :math:`Re = C|\dot{m}|`), :math:`\lambda(\dot{m})` must be an even
+    function of :math:`\dot{m}`.
+
+    Its derivative :math:`\frac{\mathrm{d}\lambda}{\mathrm{d}\dot{m}}`
+    is therefore an odd function – it changes sign when the flow direction reverses.
+
+    Since :math:`\lambda` decreases as the Reynolds number increases
+    (as seen in the Moody chart), :math:`\frac{\mathrm{d}\lambda}{\mathrm{d}\dot{m}}`
+    must be negative for positive mass flow.
+
+
+.. autoclass:: FrictionFactorModel
+    :members:
+
+
+.. _built_in_friction_factor_models:
+
+Built-in Friction Factor Models
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+There are several models provided in pandapipes by default.
+
+.. autoclass:: Nikuradse
+    :members:
+
+.. autoclass:: SwameeJain
+    :members:
+
+.. autoclass:: Colebrook
+    :members:
+
+.. autoclass:: RegimeAwareFrictionFactorModel
+    :members:
+
+
+.. _custom_friction_factor_models:
+
+Custom Friction Factor Models
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+You can implement the :class:`~FrictionFactorModel` protocol and pass it to :ref:`pipeflow function<pipeflow>`:
+
+.. code-block:: python
+
+    pp.pipeflow(net, friction_model=MyFrictionFactorModel())

src/pandapipes/__init__.py

@@ -31,5 +31,6 @@
 from pandapipes.pipeflow import *
 from pandapipes.toolbox import *
 from pandapipes.pf.pipeflow_setup import *
+from pandapipes.pf.friction_factor_model import *
 from pandapipes.std_types import *
 import pandapipes.plotting

src/pandapipes/pf/derivative_calculation.py

@@ -1,4 +1,5 @@
 import numpy as np
+
 from pandapipes.constants import NORMAL_TEMPERATURE
 from pandapipes.idx_branch import (LENGTH, D, K, RE, LAMBDA, LOAD_VEC_BRANCHES, JAC_DERIV_DM, JAC_DERIV_DP,
                                    JAC_DERIV_DP1, JAC_DERIV_DM_NODE, FROM_NODE, TO_NODE, TOUTINIT, AREA,
@@ -10,7 +11,6 @@
 from pandapipes.pf.pipeflow_setup import get_net_option, get_lookup
 from pandapipes.properties.fluids import get_fluid
 from pandapipes.properties.properties_toolbox import get_branch_real_density, get_branch_real_eta, get_branch_cp
-from scipy.optimize import newton
 
 
 def calculate_derivatives_hydraulic(net,
@@ -41,7 +41,6 @@ def calculate_derivatives_hydraulic(net,
                                                       calc_medium_pressure_with_derivative_np as calc_medium_pressure_with_derivative)
     fluid = get_fluid(net)
     gas_mode = fluid.is_gas
-    friction_model = options["friction_model"]
 
     from_nodes = branch_pit[:, FROM_NODE].astype(np.int32)
     to_nodes = branch_pit[:, TO_NODE].astype(np.int32)
@@ -57,10 +56,27 @@ def calculate_derivatives_hydraulic(net,
     eta = get_branch_real_eta(fluid, node_pit, branch_pit, p_m)
 
     # Darcy Friction factor: lambda
-    lambda_, re = calc_lambda(branch_pit[:, MDOTINIT], eta, branch_pit[:, D], branch_pit[:, K], gas_mode,
-        friction_model, branch_pit[:, LENGTH], options, branch_pit[:, AREA])
-    der_lambda = calc_der_lambda(branch_pit[:, MDOTINIT], eta, branch_pit[:, D], branch_pit[:, K], friction_model,
-                                 lambda_, branch_pit[:, AREA], re, branch_pit[:, LENGTH])
+    re = (
+        np.abs(branch_pit[:, MDOTINIT])
+        * branch_pit[:, D]
+        / (eta * branch_pit[:, AREA])
+    )
+    mask = (
+        ~np.isclose(re, 0)
+        & ~np.isclose(branch_pit[:, LENGTH], 0, rtol=1e-10, atol=1e-11)
+    )
+    k_over_D = branch_pit[mask, K] / branch_pit[mask, D]
+    lambda_ = np.zeros_like(re)
+    der_lambda = np.zeros_like(re)
+    friction_factor_model = options["friction_model"]
+    lambda_[mask], der_lambda[mask] = (
+        friction_factor_model.compute_lambda_and_dlambda_dm(
+            k_over_D,
+            re[mask],
+            branch_pit[mask, MDOTINIT],
+        )
+    )
+
     branch_pit[:, RE] = re
     branch_pit[:, LAMBDA] = lambda_
 
@@ -148,169 +164,3 @@ def get_derived_values(node_pit, from_nodes, to_nodes, use_numba):
         return calc_derived_values_numba(node_pit, from_nodes, to_nodes)
     from pandapipes.pf.derivative_toolbox import calc_derived_values_np
     return calc_derived_values_np(node_pit, from_nodes, to_nodes)
-
-
-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).
-
-    :param m:
-    :type m:
-    :param eta:
-    :type eta:
-    :param d:
-    :type d:
-    :param k:
-    :type k:
-    :param gas_mode:
-    :type gas_mode:
-    :param friction_model:
-    :type friction_model:
-    :param lengths:
-    :type lengths:
-    :param options:
-    :type options:
-    :param area:
-    :type area:
-    :return:
-    :rtype:
-    """
-    if options["use_numba"]:
-        from pandapipes.pf.derivative_toolbox_numba import (
-            calc_lambda_nikuradse_incomp_numba as calc_lambda_nikuradse_incomp,
-            calc_lambda_nikuradse_comp_numba as calc_lambda_nikuradse_comp)
-    else:
-        from pandapipes.pf.derivative_toolbox import (calc_lambda_nikuradse_incomp_np as calc_lambda_nikuradse_incomp,
-                                                      calc_lambda_nikuradse_comp_np as calc_lambda_nikuradse_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":
-        # TODO: move this import to top level if possible
-        from pandapipes.pipeflow import PipeflowNotConverged
-        max_iter = options.get("max_iter_colebrook", 100)
-        tolerance = options.get("tolerance_colebrook", 1e-4)
-        converged, lambda_colebrook = colebrook_white(re, d, k, lambda_nikuradse, max_iter, lengths, tolerance)
-        if not converged:
-            raise PipeflowNotConverged("The Colebrook-White algorithm did not converge. There might be model "
-                                       "inconsistencies. The maximum iterations can be given as 'max_iter_colebrook' "
-                                       "argument to the pipeflow.")
-        return lambda_colebrook, re
-    elif friction_model == "swamee-jain":
-        # 1.325 instead of 0.25???
-        lambda_swamee_jain = 0.25 / ((np.log10(k / (3.7 * d) + 5.74 / (re ** 0.9))) ** 2)
-        return lambda_swamee_jain, re
-    else:
-        # lambda_tot = np.where(re > 2300, lambda_laminar + lambda_nikuradse, lambda_laminar)
-        lambda_tot = lambda_laminar + lambda_nikuradse
-        return lambda_tot, re
-
-
-def calc_der_lambda(m, eta, d, k, friction_model, lambda_pipe, area, re, lengths):
-    """
-    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).
-
-    :param m:
-    :type m:
-    :param eta:
-    :type eta:
-    :param d:
-    :type d:
-    :param k:
-    :type k:
-    :param friction_model:
-    :type friction_model:
-    :param lambda_pipe:
-    :type lambda_pipe:
-    :param area:
-    :type area:
-    :return:
-    :rtype:
-    """
-
-    b_term = np.zeros_like(m)
-    df_dm = np.zeros_like(m)
-    df_dlambda = np.zeros_like(m)
-    lambda_der = np.zeros_like(m)
-    pos = ~np.isclose(re, 0)
-
-    if friction_model == "colebrook":
-        pos &= ~np.isclose(lengths, 0, rtol=1e-10, atol=1e-11)
-        b_term[pos] = (2.51 * eta[pos] * area[pos] / (m[pos] * d[pos] * np.sqrt(lambda_pipe[pos])) + k[pos] / (
-                    3.71 * d[pos]))
-
-        df_dm[pos] = -2 * 2.51 * eta[pos] * area[pos] / (m[pos] ** 2 * np.sqrt(lambda_pipe[pos]) * d[pos]) / (
-                    np.log(10) * b_term[pos])
-
-        df_dlambda[pos] = -0.5 * lambda_pipe[pos] ** (-3 / 2) - (2.51 * eta[pos] * area[pos] / (d[pos] * m[pos])) * \
-                          lambda_pipe[pos] ** (-3 / 2) / (np.log(10) * b_term[pos])
-
-        lambda_der[pos] = df_dm[pos] / df_dlambda[pos]
-
-        return lambda_der
-    elif friction_model == "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
-        # * abs(v_corr)^1.9)
-        lambda_der[pos] = 0.5 * np.log(10) ** 2 / (np.log(param) ** 3) / param * 5.166 * (
-                    (eta[pos] * area[pos]) / (d[pos])) ** 0.9 * np.abs(m[pos]) ** -1.9
-        return lambda_der
-    else:
-        lambda_der[pos] = -(64 * eta[pos] * area[pos]) / (m[pos] ** 2 * d[pos])
-        return lambda_der
-
-
-def colebrook_white(re, d, k, lambda_nikuradse, max_iter, lengths, tolerance=1e-4):
-    """
-    Function calculates the friction factor of a pipe using the Colebrook-White equation. It is an
-    implicit equation which is solved using the Newton-Raphson method. For pipes with zero flow or
-    zero length, the initial guess is returned. This should be uncritical, as the pressure loss
-    term will equal zero (lambda * u^2 * l / d).
-
-    :param re: Reynolds number [dimensionless]
-    :type re: np.array
-    :param d: Diameter [m]
-    :type d: np.array
-    :param k: Roughness [m]
-    :type k: np.array
-    :param lambda_nikuradse: Initial guess for lambda (from Nikuradse)
-    :type lambda_nikuradse: np.array
-    :param max_iter: Maximum number of iterations for the Colebrook-White calculation
-    :type max_iter: int
-    :param lengths: Length of the pipes [m] - only used to identify zero-length pipes
-    :type lengths: np.array
-    :param tolerance: Tolerance for the Colebrook-White calculation
-    :type tolerance: float
-    :return: lambda_cb, converged
-    1. lambda_cb: Friction factor according to Colebrook-White
-    2. converged: True, if the Colebrook-White calculation converged for all pipes
-    :rtype: (np.array, bool)
-    """
-
-    def colebrook_white_implicit(lambda_cb, re_nz, k_nz, d_nz):
-        return lambda_cb ** (-1 / 2) + 2 * np.log10(2.51 / (re_nz * np.sqrt(lambda_cb)) + k_nz / (3.71 * d_nz))
-
-    def cw_derivative(lambda_cb, re_nz, k_nz, d_nz):
-        return -1 / 2 * lambda_cb ** (-3 / 2) - (2.51 / re_nz) * lambda_cb ** (-3 / 2) / (
-                    np.log(10) * (2.51 / (re_nz * np.sqrt(lambda_cb)) + k_nz / (3.71 * d_nz)))
-
-    mask = ~np.isclose(re, 0) & ~np.isclose(lengths, 0, rtol=1e-10, atol=1e-11)
-    lambda_res = lambda_nikuradse
-
-    res = newton(colebrook_white_implicit, lambda_res[mask], maxiter=max_iter, args=(re[mask], k[mask], d[mask]),
-                 tol=tolerance, full_output=True, fprime=cw_derivative)  # , fprime2=cw_derivative_2)
-
-    if lambda_res[mask].size == 1:
-        lambda_res[mask] = res[0]
-        converged = res[1].converged
-    else:
-        lambda_res[mask] = res.root
-        converged = np.all(res.converged)
-
-    return converged, lambda_res