r'''Calculates Nusselt number of a fluid passing over a bed of particles
    using a correlation shown in [3]_ and cited as from [1]_ and [2]_. Likely
    the best available model as the author of [1]_ is the same as [2]_ and
    [3]_.

    .. math::
        Nu = f_a Nu_{sphere}

        Nu_{sphere} = 2 + \sqrt{Nu_{m,lam}^2 + Nu_{m,turb}^2}

        Nu_{m,lam} = 0.664Re^{0.5} Pr^{1/3}

        Nu_{m,turb} = \frac{0.037Re^{0.8} Pr}{1 + 2.443Re^{-0.1}(Pr^{2/3} -1)}

        Re = \frac{\rho v_s d_p}{\mu \epsilon}

    Parameters
    ----------
    dp : float
        Equivalent spherical particle diameter of packing [m]
    voidage : float
        Void fraction of bed packing [-]
    vs : float
        Superficial velocity of the fluid [m/s]
    rho : float
        Density of the fluid [kg/m^3]
    mu : float
        Viscosity of the fluid, [Pa*s]
    Pr : float
        Prandtl number of the fluid []
    fa : float, optional
        Fator increasing heat transfer []

    Returns
    -------
    Nu : float
        Nusselt number for heat transfer to the packed bed [-]

    Notes
    -----
    `fa` is a factor relating how much more heat transfer happens than would
    normally, around one sphere. For spheres of the same size,
    :math:`f_a = 1 + 1.5(1-\epsilon)`. For cylinders with l/d ratio of
    0.24 < l/d < 1.2 use fa = 1.6. For cubes, use fa = 1.6 For Raschig rings,
    use `fa` = 2.1 For Berl saddles, use `fa` = 2.3. fa is calculated with
    the relationship for spheres if not provided.

    Confirmed with experimental data for a range of :math:`1E-1 < Re <1,000`
    and :math:`0.4 < Pr < 1000` for spheres. Limits are smaller for other
    shapes.

    Examples
    --------
    >>> Nu_packed_bed_Gnielinski(dp=8E-4, voidage=0.4, vs=1, rho=1E3, mu=1E-3, Pr=0.7)
    61.37823202546954

    References
    ----------
    .. [1] Gnielinski, V. (1981) "Equations for the calculation of heat and
       mass transfer during flow through stationary spherical packings at
       moderate and high Peclet numbers". International Chemical Engineering
       21 (3): 378-383
    .. [2] Gnielinski, V. (1982) "Berechnung des Warmeund Stoffaustauschs in
       durchstomten ruhenden Schuttungen". Verfahrenstechnik 16(1): 36-39
    .. [3] Gnielinski, V. in G esellschaft, V. D. I., ed. VDI Heat Atlas.
       2nd ed. 2010 edition. Berlin; New York: Springer, 2010.

def Nu_packed_bed_Gnielinski(dp, voidage, vs, rho, mu, Pr, fa=None):
    r'''Calculates Nusselt number of a fluid passing over a bed of particles
    using a correlation shown in [3]_ and cited as from [1]_ and [2]_. Likely
    the best available model as the author of [1]_ is the same as [2]_ and
    [3]_.

    .. math::
        Nu = f_a Nu_{sphere}

        Nu_{sphere} = 2 + \sqrt{Nu_{m,lam}^2 + Nu_{m,turb}^2}

        Nu_{m,lam} = 0.664Re^{0.5} Pr^{1/3}

        Nu_{m,turb} = \frac{0.037Re^{0.8} Pr}{1 + 2.443Re^{-0.1}(Pr^{2/3} -1)}

        Re = \frac{\rho v_s d_p}{\mu \epsilon}

    Parameters
    ----------
    dp : float
        Equivalent spherical particle diameter of packing [m]
    voidage : float
        Void fraction of bed packing [-]
    vs : float
        Superficial velocity of the fluid [m/s]
    rho : float
        Density of the fluid [kg/m^3]
    mu : float
        Viscosity of the fluid, [Pa*s]
    Pr : float
        Prandtl number of the fluid []
    fa : float, optional
        Fator increasing heat transfer []

    Returns
    -------
    Nu : float
        Nusselt number for heat transfer to the packed bed [-]

    Notes
    -----
    `fa` is a factor relating how much more heat transfer happens than would
    normally, around one sphere. For spheres of the same size,
    :math:`f_a = 1 + 1.5(1-\epsilon)`. For cylinders with l/d ratio of
    0.24 < l/d < 1.2 use fa = 1.6. For cubes, use fa = 1.6 For Raschig rings,
    use `fa` = 2.1 For Berl saddles, use `fa` = 2.3. fa is calculated with
    the relationship for spheres if not provided.

    Confirmed with experimental data for a range of :math:`1E-1 < Re <1,000`
    and :math:`0.4 < Pr < 1000` for spheres. Limits are smaller for other
    shapes.

    Examples
    --------
    >>> Nu_packed_bed_Gnielinski(dp=8E-4, voidage=0.4, vs=1, rho=1E3, mu=1E-3, Pr=0.7)
    61.37823202546954

    References
    ----------
    .. [1] Gnielinski, V. (1981) "Equations for the calculation of heat and
       mass transfer during flow through stationary spherical packings at
       moderate and high Peclet numbers". International Chemical Engineering
       21 (3): 378-383
    .. [2] Gnielinski, V. (1982) "Berechnung des Warmeund Stoffaustauschs in
       durchstomten ruhenden Schuttungen". Verfahrenstechnik 16(1): 36-39
    .. [3] Gnielinski, V. in G esellschaft, V. D. I., ed. VDI Heat Atlas.
       2nd ed. 2010 edition. Berlin; New York: Springer, 2010.
    '''
    Re = rho*vs*dp/mu/voidage
    Nu_lam = 0.664*Re**0.5*Pr**(1/3.)
    Nu_turb = 0.037*Re**0.8*Pr/(1 + 2.443*Re**-0.1*(Pr**(2/3.)-1))
    Nu_sphere = 2 + (Nu_lam**2 + Nu_turb**2)**0.5
    if not fa:
        fa = 1 + 1.5*(1-voidage)
    return fa*Nu_sphere