v6

duct_sizer.py

"""
Duct sizing functions based on ASHRAE with extentions for Darcy/Colebrook formulas from Niazkar.
"""
# fmt: off
from datetime import date

__author__ = 'David Stein <dstein@ventrop.com>'
__date__ = date.fromisoformat('2022-11-09')
__company__ = 'Ventrop ECG, PLLC'
__all__ = ["find_hydraulic_diameter", "find_second_dimension", "get_all_possible_dims", "minimum_ventilation_rates", "run"]
# fmt: on

# region --- Imports ---------------------------------------------------------------------------------------------------
from sympy import *  # symbolical algebra library

from numbers import Real
import math

from functools import wraps
from typing import Any, TypeVar, Callable, Optional, Final

import pint
from pint import UnitRegistry  # import the unit handler

import numpy as np  # fast mathematical calculations library

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

from scipy.optimize import minimize_scalar

pd.options.display.max_columns = 20

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

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


# region --- Riser Tracking Class --------------------------------------------------------------------------------------
diversity = lambda x: np.interp(
    x,
    [0, 2, 4, 8, 10, 12, 16, 20, 24, 28, 32, 36, 40],
    [1, 1, 0.89, 0.79, 0.72, 0.64, 0.57, 0.53, 0.47, 0.44, 0.41, 0.40, 0.40],
)


class Simple_Riser:
    def __init__(self, riser_dict: dict[int, int]):
        """
        Given a dict {-1: 200, 0: 0, 1: 250, 2: 450, ...}, it assumes the airflow is given in CFM and the floor number is an int.
        """
        self.riser_dict = riser_dict
        self.branches = list(riser_dict.values())
        self.cfms = []
        for i, cfm in enumerate(self.branches):
            self.cfms.append(sum(self.branches[: i + 1]))

        self.riser_dict = {k: v for k, v in zip(riser_dict.keys(), self.cfms)}

    def run(self, dims: Optional[list[int]] = None):
        return run(self.cfms, dims)


class Diversity_Riser:
    def __init__(
        self,
        connections_dict: dict[int, int],
        cfm_per_connection: int,
        diversity_factor: Optional[float] = None,
    ):
        """
        Given a dict {-1: 2, 0: 0, 1: 3, 2: 1, ...}, it assumes the numbers given in values are number of connections and the floor number is an int.
        """
        self.connections_dict = connections_dict
        self.branches = list(connections_dict.values())
        self.total_connections = sum(self.branches)
        if diversity_factor is None:
            self.diversity_factor = diversity(sum(self.branches))
        else:
            self.diversity_factor = diversity_factor

        self.worst_case_connections = math.ceil(
            self.diversity_factor * self.total_connections
        )

        self.cfms = []
        for i, _ in enumerate(self.branches):
            if sum(self.branches[: i + 1]) <= self.worst_case_connections:
                self.cfms.append(sum(self.branches[: i + 1]) * cfm_per_connection)
            else:
                self.cfms.append(self.cfms[-1])

        self.riser_dict = {k: v for k, v in zip(connections_dict.keys(), self.cfms)}

    def run(self, dims: Optional[list[int]] = None):
        return run(self.cfms, dims)


# endregion

# region --- Custom Unit Definitions -----------------------------------------------------------------------------------
# fmt: off
ureg.define("sf = sq_ft")
ureg.define("sqft = 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)")
ureg.define("person = [people]")
ureg.define("[occupant_density] = [people] / [area]")
ureg.define("ppl_per_1000_sqft = person / (1000 * sqft)")
ureg.define("btuh = Btu / hr")
ureg.define("mbh = 1000 * btuh")
# fmt: on
# endregion ------------------------------------------------------------------------------------------------------------

# region --- Add String Support for Future GUI -------------------------------------------------------------------------
def _string_unit_input(function: F) -> Any:
    def _convert(
        x: QUANTITYLIKE | dict[str, QUANTITYLIKE]
    ) -> Q | Real | 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, Real):
            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


# endregion

# region --- Physical Constants ----------------------------------------------------------------------------------------
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")
# endregion

# region --- Defined Constants -----------------------------------------------------------------------------------------
max_vel: Final[Q] = Q("1000 fpm")
max_head_loss: Final[Q] = Q("0.1 in_WC_per_100_ft")
air_vel_tol: Final[Q] = Q("10 fpm")
# endregion

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


# endregion

# region --- Raw Equations ---------------------------------------------------------------------------------------------
@_string_unit_input
@ureg.wraps(None, ("m", None, "m"), strict=False)
def Niazkar_friction_factor(
    hydraulic_diameter: QUANTITYLIKE,
    Reynolds_number: Real,
    roughness_factor: QUANTITYLIKE = absolute_roughness,
) -> float:  # [doi:10.1007/s12205-019-2217-1]
    """
    Calculate friction factor based on pipe diameter, reynolds number, and surface conditions.
    :param hydraulic_diameter: [m]
    :param Reynolds_number: [-]
    :param roughness_factor: [m]
    :return: Niazkar friction factor [-]
    """
    if hydraulic_diameter == 0:
        return 0

    D = hydraulic_diameter
    Re = Reynolds_number
    ε = roughness_factor

    p_1 = 4.5547
    p_2 = 0.8784

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

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


@_string_unit_input
@ureg.wraps("m", ("m^3 / s", None, "Pa / m", "kg / m^3"), strict=False)
def hydraulic_diameter(
    volumetric_flow_rate: QUANTITYLIKE,
    friction_factor: float,
    head_loss: QUANTITYLIKE = max_head_loss,
    fluid_density: QUANTITYLIKE = air_density,
) -> Q:
    """
    Calculate necesary pipe diameter to achieve a volumetric flow rate based on friction factor, head loss, and fluid density.
    :param volumetric_flow_rate:  [m^3/s]
    :param friction_factor: [-]
    :param head_loss: [Pa/m]
    :param fluid_density: [kg/m^3]
    :return: Minimum pipe diameter [m]
    """
    if volumetric_flow_rate == 0:
        return 0

    f = friction_factor
    ρ = fluid_density

    return pow(
        8 * pow(volumetric_flow_rate, 2) * f * ρ / 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: QUANTITYLIKE,
    velocity: QUANTITYLIKE,
    kinematic_viscosity: QUANTITYLIKE = air_kin_viscosity,
) -> Real:
    """
    Calculate Reynold's number based on pipe diameter, fluid velocity, and fluid kinematic viscosity.
    :param hydraulic_diameter: [m]
    :param velocity: [m/s]
    :param kinematic_viscosity: [m^2/s]
    :return: Reynold's number [-]
    """
    if hydraulic_diameter == 0:
        return 0

    D_h = hydraulic_diameter
    V = velocity
    ν = kinematic_viscosity

    return D_h * V / ν


@_string_unit_input
@ureg.wraps("m/s", ("m^3 / s", "m^2"), strict=False)
def fluid_velocity(volumetric_flow_rate: QUANTITYLIKE, orifice_area: QUANTITYLIKE) -> Q:
    """
    Calculate fluid velocity through a given area based on volumetric flow rate.
    :param volumetric_flow_rate: [m^3/s]
    :param orifice_area:  [m^2]
    :return: Fluid velocity [m/s]
    """
    A = orifice_area

    return volumetric_flow_rate / A


@_string_unit_input
@ureg.wraps("m^2", "m", strict=False)
def pipe_area(internal_diameter: QUANTITYLIKE) -> Q:
    """
    Calculate pipe area based on pipe diameter.
    :param internal_diameter: [m]
    :return: Pipe area [m^2]
    """
    D = internal_diameter

    return np.pi / 4 * pow(D, 2)


@_string_unit_input
@ureg.wraps("m/s", ("m^3 / s", "m"), strict=False)
def fluid_velocity_thru_pipe(
    volumetric_flow_rate: QUANTITYLIKE, internal_diameter: QUANTITYLIKE
) -> Q:
    """
    Wrapper for :py:func:fluid_velocity().
    Calculate fluid velocity through a pipe given volumetric flow rate and pipe diameter.
    :param volumetric_flow_rate: [m^3/s]
    :param internal_diameter: [m]
    :return: Fluid velocity [m/s]
    """
    return fluid_velocity(
        volumetric_flow_rate=volumetric_flow_rate,
        orifice_area=pipe_area(internal_diameter=internal_diameter),
    )


@_string_unit_input
@ureg.wraps("inch", ("inch", "inch"), strict=False)
def equiv_diameter_rect_duct(a: QUANTITYLIKE, b: QUANTITYLIKE) -> Q:
    """
    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)


# endregion

# region --- Wrappers and User Functions -------------------------------------------------------------------------------
@_string_unit_input
@ureg.wraps("m", ("m^3/s", "Pa / m", "kg/m^3", "m", "m^2/s"), strict=False)
def find_hydraulic_diameter(
    volumetric_flow_rate: QUANTITYLIKE,
    head_loss: QUANTITYLIKE = max_head_loss,
    fluid_density: QUANTITYLIKE = air_density,
    roughness_factor: QUANTITYLIKE = absolute_roughness,
    fluid_kinematic_viscosity: QUANTITYLIKE = air_kin_viscosity,
) -> Q:
    """
    Finds hydraulic diameter based only on flow rate and headloss.

    Numerical solution of Darcy eq (verify) to determine a hydraulic diameter without a known friction factor.

    :param volumetric_flow_rate: [m^3/s]
    :param head_loss: [Pa/m]
    :param fluid_density: [kg/m^3]
    :param roughness_factor: [m]
    :param fluid_kinematic_viscosity: [m^2/s]
    :return: Hydraulic diameter [m]
    """
    if volumetric_flow_rate == 0:
        return 0

    ρ = fluid_density
    ε = roughness_factor
    ν = fluid_kinematic_viscosity

    rhs = 8 * pow(volumetric_flow_rate, 2) * ρ / pow(np.pi, 2) / head_loss

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

        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",
    ("m^3/s", "inch", "fpm", "Pa / m", "kg/m^3", "m", "m^2/s", "m/s"),
    strict=False,
)
def find_second_dimension(
    volumetric_flow_rate: QUANTITYLIKE,
    first_dimension: QUANTITYLIKE,
    max_vel: QUANTITYLIKE = max_vel,
    head_loss: QUANTITYLIKE = max_head_loss,
    fluid_density: QUANTITYLIKE = air_density,
    roughness_factor: QUANTITYLIKE = absolute_roughness,
    fluid_kinematic_viscosity: QUANTITYLIKE = air_kin_viscosity,
    abs_air_vel_tol: QUANTITYLIKE = air_vel_tol,
) -> Q:
    """
    Find second rectangular dimension of a duct based on first and flowrate/headloss.
    :param volumetric_flow_rate: [m^3/s]
    :param first_dimension: [inch]
    :param max_vel: [fpm]
    :param head_loss: [Pa/m]
    :param fluid_density: [kg/m^3]
    :param roughness_factor: [m]
    :param fluid_kinematic_viscosity: [m^2/s]
    :param abs_air_vel_tol: [m/s]
    :return: Other dimension in a rectangular duct [inch]
    """
    if volumetric_flow_rate == 0:
        return 0

    ρ = fluid_density
    ε = roughness_factor
    ν = fluid_kinematic_viscosity

    D_h = (
        find_hydraulic_diameter(
            volumetric_flow_rate=volumetric_flow_rate,
            fluid_density=ρ,
            head_loss=head_loss,
            roughness_factor=ε,
            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
    second_dimension = _round_up_to_nearest(second_dimension, 2)

    second_dimension = Q(second_dimension, "inch")

    equiv_diameter = equiv_diameter_rect_duct(first_dimension, second_dimension)

    air_vel = fluid_velocity_thru_pipe(
        volumetric_flow_rate=volumetric_flow_rate, internal_diameter=equiv_diameter
    )

    max_vel = Q(max_vel, "fpm").to("m/s")

    while (air_vel - max_vel) > Q(abs_air_vel_tol, "m/s"):
        second_dimension += Q("2 inch")
        equiv_diameter = equiv_diameter_rect_duct(first_dimension, second_dimension)

        # air_vel = fluid_velocity(
        #     volumetric_flow_rate=flow_rate, orifice_area=first_dim * second_dim
        # )
        air_vel = fluid_velocity_thru_pipe(
            volumetric_flow_rate=volumetric_flow_rate, internal_diameter=equiv_diameter
        )

    return second_dimension


@_string_unit_input
@ureg.wraps(None, ("cfm", "inch", None, None), strict=False)
def get_all_possible_dims(
    volumetric_flow_rate: QUANTITYLIKE,
    max_dim: QUANTITYLIKE = Q("40 inch"),
    max_aspect_ratio: Real = 3.5,
    print_: bool = True,
):
    Q_dot = Q(volumetric_flow_rate, "cfm")
    dims = []
    for i in [5] + list(range(6, max_dim + 1))[::2]:
        second_dimension = find_second_dimension(Q_dot, i)

        if second_dimension.m >= i:
            AR = second_dimension.m / i
        else:
            AR = i / second_dimension.m

        if AR <= max_aspect_ratio:

            dims.append((Q(i, "inch"), second_dimension))

            if print_:
                if AR == 1:
                    print(i, "x", second_dimension.m, " *")
                else:
                    print(i, "x", second_dimension.m)

    return dims


@_string_unit_input
@ureg.wraps(
    ("cfm", "cfm", "person"),
    (
        "sqft",
        "cfm / person",
        "cfm / sqft",
        "ppl_per_1000_sqft",
        "cfm / sqft",
        "person",
        "sqft / person",
    ),
    strict=False,
)
def minimum_ventilation_rates(
    zone_floor_area: QUANTITYLIKE,
    people_outdoor_air_rate: QUANTITYLIKE,
    area_outdoor_air_rate: QUANTITYLIKE,
    default_occupant_density: QUANTITYLIKE,
    exhaust_airflow_rate: QUANTITYLIKE,
    exact_occupancy: QUANTITYLIKE = Q("0 person"),
    exact_area_per_person: QUANTITYLIKE = Q("0 sqft / person"),
):
    A_z = zone_floor_area
    R_p = people_outdoor_air_rate
    R_a = area_outdoor_air_rate

    if exact_occupancy:
        P_z = math.ceil(exact_occupancy)
    elif exact_area_per_person:
        P_z = math.ceil(A_z / exact_area_per_person)
    elif default_occupant_density:
        P_z = math.ceil(A_z * default_occupant_density / 1000)
    else:
        P_z = 0

    V_bz = _round_up_to_nearest(R_p * P_z + R_a * A_z, 10)

    exhaust_airflow = exhaust_airflow_rate * A_z

    if exhaust_airflow == 0:
        print(f"{V_bz} cfm OA/RA, {P_z} people")
    else:
        exhaust_airflow = _round_up_to_nearest(exhaust_airflow, 10)
        print(f"{V_bz} cfm OA, {exhaust_airflow} cfm EX, {P_z} people")

    return V_bz, exhaust_airflow, P_z


# endregion


def run(cfms: Optional[list[int]] = None, first_dims: Optional[list[int]] = None):
    if cfms is None:
        cfms = [int(x) for x in input("cfms: ").split()]

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

    df["D_h [inch]"] = df["Flow [cfm]"].apply(
        lambda x: find_hydraulic_diameter(Q(x, "cfm")).to("inch").m
    )

    if first_dims is None:
        print(df)
        first_dims = [int(x) for x in input("first dims: ").split()]

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

        df["D_e [inch]"] = df.apply(
            lambda x: equiv_diameter_rect_duct(
                Q(x["Width [inch]"], "inch"), Q(x["Height [inch]"], "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_thru_pipe(
                Q(x["Flow [cfm]"], "cfm"), Q(x["D_e [inch]"], "inch")
            )
            .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")
            ),
            axis=1,
        )
        df["Friction Factor"] = df.apply(
            lambda x: Niazkar_friction_factor(
                Q(x["D_e [inch]"], "inch"), x["Reynolds Number"]
            ),
            axis=1,
        )
        df["Velocity pressure [in WC]"] = df.apply(
            lambda x: (air_density * pow(Q(x["Fluid velocity [fpm]"], "fpm"), 2) / 2)
            .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(math.ceil)
        df["Reynolds Number"] = df["Reynolds Number"].apply(math.ceil)
        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)
    return df

    # while True:
    #     first_dims = [int(x) for x in input("first dims: ").split()]
    #
    #     calc(df, first_dims)
    #
    #     print(df)