Implemented new friction factor calculation spanning entire flow regime

[cvTHM]

The changes include the implementation of a new calculation method for the pipe friction factor for Darcy-Weisbach equation. It spans the entire flow regime for laminar and turbulent flows and the transition region and is an explicit expression.
Total deviations from the Swamee-Jain equation in the turbulent regime and from the solely laminar approach 64/Re can be found in the attached image.

The original publication for this correlation is:
Churchill, S. W. Friction-factor equation spans all fluid-flow regimes, 1977, Chemical Engineering, p. 91-92

[scratch_test_churchill_friction_factor.txt](https://github.com/user-attachments/files/23544694/scratch_test_churchill_friction_factor.txt)

[cvTHM]

Might be related to a recent question raised in #754

[EPrade]

Hey looks good. Could you maybe add a test case to validate the results with literature and a changelog entry?

[cvTHM]

Hi,
I added an entry to the changelog file. Where in the file directory should I put the test case? The text file attached to the initial commit contains the code snippet below (comparison to Swamee-Jain friction factor and creation of the plot above). Should It put this into a separate py file?

# %%
import numpy as np
from scipy.optimize import newton


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":
        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 == "churchill":
        paramA = (-2.457*np.log((7/re)**0.9 + 0.27*k/d))**16
        paramB = (37530/re)**16
        lambda_churchill = 8*((8/re)**12 + 1/(paramA+paramB)**1.5)**(1/12)
        return lambda_churchill, re

    else:
        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
    
    elif friction_model == "churchill":
        param = (7*eta[pos]*area[pos]/(m[pos]*d[pos]))**0.9 + 0.27*k[pos]/d[pos]
        partial_dparamdm = -0.9*7**0.9 * (d[pos]/(eta[pos] * area[pos]))**(-0.9) * m[pos]**(-1.9)

        paramsAB = ((-2.457*np.log((7*eta[pos]*area[pos]/(m[pos] * d[pos]))**0.9 + 0.27*k[pos]/d[pos]))**16 + (37530*eta[pos]*area[pos]/(m[pos] * d[pos]))**16)
        paramC = (8*eta[pos]*area[pos]/(m[pos]*d[pos]))**12 + paramsAB**(-1.5)
        
        partial_dAdm = 16*(2.457*np.log(param))**15 * 2.457*param**(-1) * partial_dparamdm
        partial_dBdm = -16*37530**16*(eta[pos]*area[pos]/(m[pos]*d[pos]))**17

        partial_dCdm = -12*(8*eta[pos]*area[pos]/(m[pos]*d[pos]))**11 * (8*eta[pos]*area[pos]/(m[pos]**2*d[pos])) - paramsAB**(-3)*1.5*paramsAB**0.5 * (partial_dAdm + partial_dBdm)

        lambda_der[pos] = 2/3* paramC**(-11/12) * partial_dCdm

        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


if __name__ == '__main__':

    import matplotlib.pyplot as plt

    re = np.linspace(1e02, 1e08, 1000000)# m * d / (eta * area)
    re_lam = np.where(re<2300)
    re_turb = np.where(re>5e03)

    d = 0.25 * np.ones(len(re))
    k = 0.0469*1e-03 * np.ones(len(re))
    eta = 1e-03 * np.ones(len(re))
    area = d**2 / 4 * np.pi
    m = re * eta * area / d

    lengths = 1 * np.ones(len(re))
    options = {'use_numba':False}

    lambda_pipe = np.ones(len(re))
   

    # Plotting of Churchill friction factor for smooth and rough pipes
    varepsilon = np.linspace(0.0000002, 0.001, 3) # Roughness (m)

    f_friction = {}
    for vareps in varepsilon:
        lam_churchill, _ = calc_lambda(m = m, eta = eta, d = d, k = vareps, area = area, lengths = lengths, options = options, gas_mode=False, friction_model = 'churchill')     

        f_friction[vareps] = lam_churchill


    fig1, ax1 = plt.subplots(layout = "tight")

    ax1.loglog(re, f_friction[varepsilon[0]], color="#007732", label='smooth')
    ax1.loglog(re, f_friction[varepsilon[-1]], color="#994D06", label='rough')
    ax1.axvline(x=2300, color='r', linestyle='--')    
    ax1.set_xlim(5e02, 1e08)
    ax1.set_ylim(5e-03, 2e-01)
    ax1.grid(True, which='both')
    ax1.set_xlabel('Reynolds number')
    ax1.set_ylabel('Friction factor')
    ax1.set_title('Churchill friction factor for smooth and rough pipes')
    ax1.legend()
    plt.show()


    # Plotting of relative deviations to other friction factors
    friction_lam = np.array([64/Re for Re in re[re_lam]]).flatten()

    lam_churchill_smooth, _ = calc_lambda(m = m, eta = eta, d = d, k = varepsilon[0], area = area, lengths = lengths, options = options, gas_mode=False, friction_model = 'churchill')

    lam_churchill_rough, _ =  calc_lambda(m = m, eta = eta, d = d, k = varepsilon[-1], area = area, lengths = lengths, options = options, gas_mode=False, friction_model = 'churchill')

    lam_swamee_jain_smooth, _ = calc_lambda(m = m, eta = eta, d = d, k = varepsilon[0], area = area, lengths = lengths, options = options, gas_mode=False, friction_model = 'swamee-jain')

    lam_swamee_jain_rough, _ = calc_lambda(m = m, eta = eta, d = d, k = varepsilon[-1], area = area, lengths = lengths, options = options, gas_mode=False, friction_model = 'swamee-jain')


    fig2, ax2 = plt.subplots(layout = "tight")

    ax2.plot(re[re_lam], ((lam_churchill_rough[re_lam]-friction_lam)/friction_lam) * 100, color="#994D06")

    ax2.plot(re[re_lam], ((lam_churchill_smooth[re_lam]-friction_lam)/friction_lam) * 100, color="#007732")

    ax2.plot(re[re_turb], ((lam_churchill_rough[re_turb]-lam_swamee_jain_rough[re_turb])/lam_swamee_jain_rough[re_turb]) * 100, color="#994D06", label=r"${\Delta}$f(rough)")

    ax2.plot(re[re_turb], ((lam_churchill_smooth[re_turb]-lam_swamee_jain_smooth[re_turb])/lam_swamee_jain_smooth[re_turb]) * 100, color="#007732", label=r"${\Delta}$f(smooth)")

    ax2.grid(True, which='both')
    ax2.set_xscale("log")
    ax2.set_xlabel('Reynolds number [-]')
    ax2.set_ylabel(f'Deviation of Churchill friction factor\nto Swamee-Jain and 64/Re [%]')
    ax2.set_ylim([-0.06, 0.15])
    ax2.legend(loc = "best")
    plt.show()