For sizing HVAC rectangular ducts (assumes galvanized steel)

duct_sizer.py

from sympy import *  # symbolical algebra library

from enum import Enum, auto
from typing import *
from numbers import Number
from fractions import Fraction
import math

import inspect
from functools import partial, wraps
from typing import Union, Any, TypeVar, Callable, cast, overload, Optional

F = TypeVar("F", bound=Callable)


import pint
from pint import UnitRegistry  # import the unit handler

ureg = UnitRegistry(
    auto_reduce_dimensions=True, non_int_type=Fraction
)  # initialize a registry to hold the unit conversions
Q = ureg.Quantity  # alias `Q` to the quantity constructor
QUANTITYLIKE = Q | str | Number

import numpy as np  # fast mathematical calculations library

from CoolProp.CoolProp import PropsSI  # fluid properties charts

import pandas as pd  # two-dimensional array (spreadsheets) library

from scipy.optimize import root, minimize_scalar

from functools import partial

pd.options.display.max_columns = 20

ureg = UnitRegistry()
Q = ureg.Quantity

ureg.define("sf = sq_ft")
ureg.define("cf = ft^3")
ureg.define("cfm = ft^3 / min")
ureg.define("one = []")
ureg.define(
    "percent = one / 100"
)  # needs to define percent in the prompts as "9 percent", cannot use symbol
ureg.define("cfs = ft^3 / second")
ureg.define("fpm = ft / min")

ureg.define("in_WC = inch_H2O_39F")
ureg.define("in_WC_per_100_ft = inch_H2O_39F / (100 * ft)")


def _string_unit_input(
    function: F,
) -> Any:
    def _convert(
        x: QUANTITYLIKE | dict[str, QUANTITYLIKE]
    ) -> Q | Number | dict[str, QUANTITYLIKE]:
        if isinstance(x, str):
            try:
                return Q(x)
            except pint.errors.UndefinedUnitError:
                raise ValueError(f"Can't convert '{x}' to a quantity.")
        elif isinstance(x, Number):
            return x
        elif isinstance(x, Q):
            return x
        elif isinstance(x, dict):
            for k, v in x.items():
                x[k] = _convert(v)
            return x

    @wraps(function)
    def _wrapper(*args: Any, **kwargs: Any) -> Any:
        args = [_convert(x) for x in args]
        kwargs = _convert(kwargs)

        return function(*args, **kwargs)

    return _wrapper


air_density: Final[Q] = Q("0.075 lb / ft^3")
_air_viscosity: Final[Q] = Q("0.0432 lb / ft / hr")
air_kin_viscosity: Final[Q] = _air_viscosity / air_density
specific_heat: Final[Q] = Q("0.24 Btu / lb / delta_degF")
energy_factor: Final[Q] = Q("1.08 Btu / hr / delta_degF / cfm")
absolute_roughness: Final[Q] = Q("0.0003 ft")
g_c: Final[Q] = Q(32.174, "(lb * ft) / (lbf * s^2)")

D_h: Final[Symbol] = symbols("D_h", positive=True)
f: Final[Symbol] = symbols("f", positive=True)
rho: Final[Symbol] = symbols("rho", positive=True)
Re: Final[Symbol] = symbols(r"\text{Re}", positive=True)
epsilon: Final[Symbol] = symbols(r"\varepsilon", positive=True)
HL: Final[Symbol] = symbols("H_L", positive=True)
V: Final[Symbol] = symbols("V", positive=True)
flow: Final[Symbol] = symbols(r"\text{(flow)}", positive=True)
A: Final[Symbol] = symbols("A", positive=True)
mu: Final[Symbol] = symbols(r"\mu", positive=True)

area = Eq(A, pi / 4 * pow(D_h, 2))
velocity = Eq(V, flow / area.rhs)
darcy = Eq(HL, f * rho * pow(velocity.rhs, 2) / (2 * D_h))


# calc_hydraulic_diameter = lambdify([HL, rho, f, flow], solve(darcy, D_h)[0], 'numpy')


# def calc_pipe_area(hydraulic_diameter: Q) -> Q:  # [area]
#     if not hydraulic_diameter.check("[length]"):
#         raise ValueError
#     return np.pi / 4 * pow(hydraulic_diameter, 2)
#
#
# def calc_velocity(
#     volumetric_flow_rate: Q, hydraulic_diameter: Optional[Q], pipe_area: Optional[Q]
# ) -> Q:  # [velocity]
#     if not volumetric_flow_rate.check("[volume] / [time]"):
#         raise ValueError
#
#     if (hydraulic_diameter is None) and (pipe_area is None):
#         raise ValueError
#     elif hydraulic_diameter is not None:
#         if not hydraulic_diameter.check("[length]"):
#             raise ValueError
#         return volumetric_flow_rate / calc_pipe_area(hydraulic_diameter)
#     elif pipe_area is not None:
#         if not pipe_area.check("[area]"):
#             raise ValueError
#         return volumetric_flow_rate / pipe_area
#
#
# def calc_hydraulic_diameter(
#     head_loss: Q, fluid_density: Q, volumetric_flow_rate: Q, friction_factor: Number
# ):
#     if not all(
#         [
#             head_loss.check("[pressure] / [length]"),
#             fluid_density.check("[density]"),
#             volumetric_flow_rate.check("[volume] / [time]"),
#         ]
#     ):
#         raise ValueError
#     return Q(
#         pow(
#             8
#             * pow(volumetric_flow_rate, 2)
#             * friction_factor
#             * fluid_density
#             / pow(np.pi, 2)
#             / head_loss,
#             1 / 5,
#         )
#         .to_base_units()
#         .m,
#         "m",
#     )
#
#
# def calc_reynolds_number(
#     head_loss: Q,
#     fluid_density: Q,
#     volumetric_flow_rate: Q,
#     friction_factor: Number,
#     kinematic_viscosity: Q,
# ) -> Number:  # dimensionless
#     if not all(
#         [
#             head_loss.check("[pressure] / [length]"),
#             fluid_density.check("[density]"),
#             volumetric_flow_rate.check("[volume] / [time]"),
#             kinematic_viscosity.check("[kinematic_viscosity]"),
#         ]
#     ):
#         raise ValueError
#     D_h = hydraulic_diameter(
#         head_loss, fluid_density, volumetric_flow_rate, friction_factor
#     )
#     V = calc_velocity(volumetric_flow_rate, hydraulic_diameter=D_h, pipe_area=None)
#
#     return float(D_h * V / kinematic_viscosity)
#
#
# def calc_friction_factor(
#     head_loss: Q,
#     fluid_density: Q,
#     volumetric_flow_rate: Q,
#     material_roughness_factor: Q,
#     kinematic_viscosity: Q,
# ) -> Number:
#     if not all(
#         [
#             head_loss.check("[pressure] / [length]"),
#             fluid_density.check("[density]"),
#             volumetric_flow_rate.check("[volume] / [time]"),
#             material_roughness_factor.check("[length]"),
#             kinematic_viscosity.check("[kinematic_viscosity]"),
#         ]
#     ):
#         raise ValueError
#
#     def f(x):
#         D_h = hydraulic_diameter(
#             head_loss=head_loss,
#             fluid_density=fluid_density,
#             volumetric_flow_rate=volumetric_flow_rate,
#             friction_factor=x,
#         )
#         Re = calc_reynolds_number(
#             head_loss=head_loss,
#             fluid_density=fluid_density,
#             volumetric_flow_rate=volumetric_flow_rate,
#             friction_factor=x,
#             kinematic_viscosity=kinematic_viscosity,
#         )
#
#             tmp_A = float(material_roughness_factor / (3.7 * D_h))
#             tmp_B = 2.51 / (Re * np.sqrt(x))
#
#             return -2 * np.log10(tmp_A + tmp_B) - 1 / np.sqrt(x)
#
#         return root(f, 0.02).x[0]
#
#
#     def get_friction_factor_simple(
#         volumetric_flow_rate: Q, head_loss: Q = Q("0.1 in_WC_per_100_ft")
#     ) -> Number:
#         return calc_friction_factor(
#             head_loss=head_loss,
#             fluid_density=air_density,
#             volumetric_flow_rate=volumetric_flow_rate,
#             material_roughness_factor=absolute_roughness,
#             kinematic_viscosity=air_kin_viscosity,
#         )
#
#
#     def get_hydraulic_diameter_simple(
#         volumetric_flow_rate: Q, head_loss: Q = Q("0.1 in_WC_per_100_ft")
#     ) -> Q:
#         friction_factor = get_friction_factor_simple(
#             volumetric_flow_rate=volumetric_flow_rate, head_loss=head_loss
#         )
#         return hydraulic_diameter(
#             head_loss=head_loss,
#             fluid_density=air_density,
#             volumetric_flow_rate=volumetric_flow_rate,
#             friction_factor=friction_factor,
#         )
#
#
#     # reynolds = Eq(
#     #     Re, solve(darcy, D_h)[0] * velocity.rhs.subs(D_h, solve(darcy, D_h)[0]) / mu
#     # )
#
#     # colebrook = Eq(
#     #     1 / sqrt(f),
#     #     -2 * log(epsilon / (3.7 * solve(darcy, D_h)[0]) + 2.51 / (Re * sqrt(f)), 10),
#     # ).subs(Re, reynolds.rhs)
#
#
#     # colebrook_func = lambdify([f, epsilon, rho, HL, Re], colebrook.rhs - colebrook.lhs)
#
#
#     def colebrook_func(f: Number, head_loss: Q, epsilon: Q, flow: Q, mu: Q, rho: Q):
#         head_loss = head_loss.to_base_units().m
#         epsilon = epsilon.to_base_units().m
#         flow = flow.to_base_units().m
#         mu = mu.to_base_units().m
#         rho = rho.to_base_units().m
#         return -2 * np.log10(
#             (5 / 37)
#             * pow(2 * np.pi, 2 / 5)
#             * pow(head_loss, 1 / 5)
#             * epsilon
#             / (pow(flow, 2 / 5) * pow(f * rho, 1 / 5))
#             + 32
#             / 51
#             * pow(2 * np.pi, 3 / 5)
#             * mu
#             * pow(rho, 1 / 5)
#             / (pow(head_loss, 1 / 5) * pow(flow, 3 / 5) * pow(f, 3 / 10))
#         ) - pow(f, -1 / 2)


# print root(f, 0.02)


@_string_unit_input
@ureg.wraps(None, ("m", None, "m"), strict=False)
def Niazkar_friction_factor(
    hydraulic_diameter: Number, Reynolds_number: Number, roughness_factor: Number
) -> Number:
    """
    Assumes the hydraulic_diameter and roughness_factor are given in meters.
    :param hydraulic_diameter: [m]
    :param Reynolds_number: [-]
    :param roughness_factor: [m]
    :return: [-]
    """

    D = hydraulic_diameter
    Re = Reynolds_number
    epsilon = roughness_factor

    p_1 = 4.5547
    p_2 = 0.8784

    A = -2 * np.log10(epsilon / (3.7 * D) + p_1 / pow(Re, p_2))
    B = -2 * np.log10(epsilon / (3.7 * D) + 2.51 * A / Re)
    C = -2 * np.log10(epsilon / (3.7 * D) + 2.51 * B / Re)

    f = pow(A - pow(B - A, 2) / (C - 2 * B + A), -2)
    return f


@_string_unit_input
@ureg.wraps("m", ("Pa / m", "kg / m^3", "m^3 / s", None), strict=False)
def hydraulic_diameter(
    head_loss: Number,
    fluid_density: Number,
    volumetric_flow_rate: Number,
    friction_factor: Number,
) -> Number:
    """
    Assumes head_loss is given in Pascals per meter, fluid_density is given in kilogram / cubic meter,
    volumetric_flow_rate is given in meters cubed per second, and friction_factor is dimensionless.
    :param head_loss:  [Pa/m]
    :param fluid_density: [kg/m^3]
    :param volumetric_flow_rate: [m^3/s]
    :param friction_factor: [-]
    :return: [m]
    """
    return pow(
        8
        * pow(volumetric_flow_rate, 2)
        * friction_factor
        * fluid_density
        / pow(np.pi, 2)
        / head_loss,
        1 / 5,
    )


@_string_unit_input
@ureg.wraps(None, ("m", "m/s", "m^2 / s"), strict=False)
def Reynolds_number(
    hydraulic_diameter: Number, velocity: Number, kinematic_viscosity: Number
) -> Number:
    """
    Assumes hydraulic_diameter given in meters, velocity in meters per second, and kinematic viscosity in meters squared per second.
    :param hydraulic_diameter: [m]
    :param velocity: [m/s]
    :param kinematic_viscosity: [m^2/s]
    :return:
    """
    D_h = hydraulic_diameter
    V = velocity
    nu = kinematic_viscosity

    return D_h * V / nu


@_string_unit_input
@ureg.wraps("m/s", ("m^3 / s", "m^2"), strict=False)
def fluid_velocity(volumetric_flow_rate: Number, orifice_area: Number) -> Number:
    """
    Assumes volumetric_flow_rate is given in cubic meters per second, and orifice_area is given in square meters.
    :param volumetric_flow_rate: [m^3 / s]
    :param orifice_area:  [m^2]
    :return: [m/s]
    """
    return volumetric_flow_rate / orifice_area


@_string_unit_input
@ureg.wraps("m^2", "m", strict=False)
def pipe_area(internal_diameter: Number) -> Number:
    """
    Assumes internal_diameter is given in meters.
    :param internal_diameter: [m]
    :return: [m^2]
    """
    return np.pi / 4 * pow(internal_diameter, 2)


@_string_unit_input
@ureg.wraps(None, ("m^3 / s", "m"), strict=False)
def fluid_velocity_thru_pipe(
    volumetric_flow_rate: Number, internal_diameter: Number
) -> Number:
    return fluid_velocity(
        volumetric_flow_rate=volumetric_flow_rate,
        orifice_area=pipe_area(internal_diameter=internal_diameter),
    )


@_string_unit_input
@ureg.wraps("m", ("m^3/s", "kg/m^3", "Pa/m", "m", "m^2/s"), strict=False)
def find_hydraulic_diameter(
    volumetric_flow_rate: Number,
    fluid_density: Number,
    head_loss: Number,
    roughness_factor: Number,
    fluid_kinematic_viscosity: Number,
) -> Number:
    """
    Finds hydraulic diameter based only on flow rate and headloss.
    :param volumetric_flow_rate: [m^3/s]
    :param fluid_density: [kg/m^3]
    :param head_loss: [Pa/m]
    :param roughness_factor: [m]
    :param fluid_kinematic_viscosity: [m^2/s]
    :return: [m]
    """
    F = volumetric_flow_rate
    rho = fluid_density
    H_L = head_loss
    epsilon = roughness_factor
    nu = fluid_kinematic_viscosity

    rhs = 8 * pow(F, 2) * rho / pow(np.pi, 2) / H_L

    def f(x):
        V = fluid_velocity_thru_pipe(volumetric_flow_rate=F, internal_diameter=x)
        Re = Reynolds_number(hydraulic_diameter=x, velocity=V, kinematic_viscosity=nu)
        N_f = Niazkar_friction_factor(
            hydraulic_diameter=x, Reynolds_number=Re, roughness_factor=epsilon
        )

        return abs(rhs - pow(x, 5) / N_f)

    return minimize_scalar(f, method="bounded", bounds=(0.001, 3)).x


@_string_unit_input
@ureg.wraps("inch", ("inch", "inch"), strict=False)
def equiv_diameter_rect_duct(a: Number, b: Number) -> Number:
    """
    Finds equivalent round duct given length and width of a rectangular duct (all in inches)
    :param a: [inch]
    :param b: [inch]
    :return: [inch]
    """
    return 1.30 * pow(a * b, 0.625) / pow(a + b, 0.250)


@_string_unit_input
@ureg.wraps("inch", ("m^3/s", "kg/m^3", "Pa/m", "m", "m^2/s", "inch"), strict=False)
def find_second_dimension(
    volumetric_flow_rate: Number,
    fluid_density: Number,
    head_loss: Number,
    roughness_factor: Number,
    fluid_kinematic_viscosity: Number,
    first_dimension: Number,
) -> Number:
    """
    Find second rectangular dimension of a duct based on first and flowrate/headloss.
    :param volumetric_flow_rate: [m^3/s]
    :param fluid_density: [kg/m^3]
    :param head_loss: [Pa/m]
    :param roughness_factor: [m]
    :param fluid_kinematic_viscosity: [m^2/s]
    :param first_dimension: [inch]
    :return: [inch]
    """
    D_h = (
        find_hydraulic_diameter(
            volumetric_flow_rate=volumetric_flow_rate,
            fluid_density=fluid_density,
            head_loss=head_loss,
            roughness_factor=roughness_factor,
            fluid_kinematic_viscosity=fluid_kinematic_viscosity,
        )
        .to("inch")
        .m
    )

    def f(x):
        return abs(D_h - equiv_diameter_rect_duct(a=first_dimension, b=x).m)

    second_dimension = minimize_scalar(f, method="bounded", bounds=(1, 200)).x

    return _round_up_to_nearest(second_dimension, 2)


def _round_up_to_nearest(x: float | int, multiple: int) -> int:
    if x < multiple:
        return multiple
    elif math.isclose(x % multiple, 0, abs_tol=0.1):
        return math.floor(x)
    return math.ceil(x / multiple) * multiple


@ureg.wraps("inch", "cfm", strict=False)
def get_hydraulic_diameter(flow_rate) -> Q:
    """

    :param flow_rate: [cfm]
    :return: [inch]
    """
    flow_rate = Q(flow_rate, "cfm")
    return find_hydraulic_diameter(
        volumetric_flow_rate=flow_rate,
        fluid_density=air_density,
        head_loss="0.1 in_WC_per_100_ft",
        roughness_factor=absolute_roughness,
        fluid_kinematic_viscosity=air_kin_viscosity,
    )


@ureg.wraps("inch", ("m^3/s", "inch", "m/s", "Pa/m", "m", "m^2/s"), strict=False)
def get_second_dim(
    flow_rate,
    first_dim,
    max_vel,
    head_loss,
    roughness_factor,
    fluid_kinematic_viscosity,
) -> Q:
    """

    :param flow_rate: [m/s]
    :param first_dim: [inch]
    :param max_vel: [m/s]
    :param head_loss: [Pa/m]
    :param roughness_factor: [m]
    :param fluid_kinematic_viscosity: [m^2/s]
    :return: [inch]
    """
    first_dim = Q(first_dim, "inch")
    max_vel = Q(max_vel, "m/s")
    second_dim = find_second_dimension(
        volumetric_flow_rate=flow_rate,
        fluid_density=air_density,
        head_loss=head_loss,
        roughness_factor=roughness_factor,
        fluid_kinematic_viscosity=fluid_kinematic_viscosity,
        first_dimension=first_dim,
    )

    air_vel = fluid_velocity(
        volumetric_flow_rate=flow_rate, orifice_area=first_dim * second_dim
    )

    while air_vel > max_vel:
        second_dim += Q("2 inch")
        air_vel = fluid_velocity(
            volumetric_flow_rate=flow_rate, orifice_area=first_dim * second_dim
        )

    return second_dim


@ureg.wraps("inch", ("cfm", "inch", "fpm", "in_WC_per_100_ft"), strict=False)
def get_second_dim_simple(flow_rate, first_dim, max_vel=None, head_loss=None):
    """

    :param flow_rate: [cfm]
    :param first_dim: [inch]
    :param max_vel: Optional[fpm]
    :param head_loss: Optional[in_WC_per_100_ft]
    :return:
    """
    if max_vel is None:
        max_vel = Q("1200 fpm")
    else:
        max_vel = Q(max_vel, "fpm")

    if head_loss is None:
        head_loss = Q("0.1 in_WC_per_100_ft")
    else:
        head_loss = Q(head_loss, "in_WC_per_100_ft")

    flow_rate = Q(flow_rate, "cfm")
    first_dim = Q(first_dim, "inch")

    return get_second_dim(
        flow_rate, first_dim, max_vel, head_loss, absolute_roughness, air_kin_viscosity
    )


if __name__ == "__main__":
    cfm = input("cfms: ")
    cfms = [int(x) for x in cfm.split()]

    df = pd.DataFrame({"Flow [cfm]": cfms})

    df["D_h [inch]"] = df["Flow [cfm]"].apply(lambda x: get_hydraulic_diameter(x).m)

    print(df)

    first_dim = input("first dims: ")
    first_dims = [int(x) for x in first_dim.split()]

    def calc(df, first_dims):
        df["Width [inch]"] = first_dims
        df["Height [inch]"] = df.apply(
            lambda x: get_second_dim_simple(x["Flow [cfm]"], x["Width [inch]"]).m,
            axis=1,
        )

        df["D_e [inch]"] = df.apply(
            lambda x: equiv_diameter_rect_duct(x["Width [inch]"], x["Height [inch]"]).m,
            axis=1,
        )
        df["Flow area [ft^2]"] = df.apply(
            lambda x: (Q(x["Width [inch]"], "inch") * Q(x["Height [inch]"], "inch"))
            .to("ft^2")
            .m,
            axis=1,
        )
        df["Fluid velocity [fpm]"] = df.apply(
            lambda x: fluid_velocity(
                Q(x["Flow [cfm]"], "cfm"), Q(x["Flow area [ft^2]"], "ft^2")
            )
            .to("fpm")
            .m,
            axis=1,
        )
        df["Reynolds Number"] = df.apply(
            lambda x: Reynolds_number(
                Q(x["D_e [inch]"], "inch"),
                Q(x["Fluid velocity [fpm]"], "fpm"),
                air_kin_viscosity,
            ),
            axis=1,
        )
        df["Friction Factor"] = df.apply(
            lambda x: Niazkar_friction_factor(
                Q(x["D_h [inch]"], "inch"), x["Reynolds Number"], absolute_roughness
            ),
            axis=1,
        )
        df["Velocity pressure [in WC]"] = df.apply(
            lambda x: (
                air_density * pow(Q(x["Fluid velocity [fpm]"], "fpm"), 2) / (2 * g_c)
            )
            .to("in_WC")
            .m,
            axis=1,
        )
        df["Head loss [in WC / 100 ft]"] = df.apply(
            lambda x: (
                pow(Q(x["Fluid velocity [fpm]"], "fpm"), 2)
                * x["Friction Factor"]
                * air_density
                / 2
                / Q(x["D_e [inch]"], "inch")
            )
            .to("in_WC_per_100_ft")
            .m,
            axis=1,
        )

        df["D_e [inch]"] = round(df["D_e [inch]"], 3)
        df["Flow area [ft^2]"] = round(df["Flow area [ft^2]"], 4)
        df["Fluid velocity [fpm]"] = df["Fluid velocity [fpm]"].apply(
            lambda x: math.ceil(x)
        )
        df["Reynolds Number"] = df["Reynolds Number"].apply(lambda x: math.ceil(x))
        df["Friction Factor"] = round(df["Friction Factor"], 5)
        df["Velocity pressure [in WC]"] = round(df["Velocity pressure [in WC]"], 4)
        df["Head loss [in WC / 100 ft]"] = round(df["Head loss [in WC / 100 ft]"], 4)

    calc(df, first_dims)
    print(df)

    while True:
        first_dim = input("first dims: ")
        first_dims = [int(x) for x in first_dim.split()]

        calc(df, first_dims)

        print(df)

Input is a list of airflows (in cfm, just type the numbers with spaces in between) from cellar floor up to roof. Second input is a guess at one side of the ducts based on hydraulic dia (in inches, just type the numbers with spaces in between).