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 fractions import Fraction

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.88, 0.78, 0.7, 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(self.riser_dict.values())
        self.cfms = []
        for i, cfm in enumerate(self.branches):
            self.cfms.append(sum(self.branches[: i + 1]))


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

    def iupdate(self, riser_dict: dict[int, int]):
        self.riser_dict |= riser_dict
        self.branches = list(self.riser_dict.values())
        self.cfms = []
        for i, cfm in enumerate(self.branches):
            self.cfms.append(sum(self.branches[: i + 1]))

    def update(self, riser_dict: dict[int, int]):
        self.iupdate(riser_dict)
        return self

    def __add__(self, other: "Simple_Riser"):
        new_riser_dict = dict()
        for floor_number in self.riser_dict:
            if floor_number in other.riser_dict:
                new_riser_dict |= {floor_number: self.riser_dict[floor_number] + other.riser_dict[floor_number]}

        return Simple_Riser(new_riser_dict)


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: list[int] = []
        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.worst_case_connections * cfm_per_connection)

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

    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("1200 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_int(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

def _round_up_to_nearest_float(x: float | int, factor: Fraction) -> float:
    if x < factor:
        return factor
    else:
        return float(math.ceil(x / float(factor)) * factor)

# 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_int(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", "fpm", None, None), strict=False)
def get_all_possible_dims(
    volumetric_flow_rate: QUANTITYLIKE,
    max_dim: QUANTITYLIKE = Q("40 inch"),
    max_vel: QUANTITYLIKE = max_vel,
    max_aspect_ratio: Real = 3.5,
    print_: bool = True,
):
    Q_dot = Q(volumetric_flow_rate, "cfm")
    dims = []
    used = []
    for i in [5] + list(range(6, max_dim + 1))[::2]:
        second_dimension = find_second_dimension(Q_dot, i, max_vel=max_vel)

        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 (second_dimension.m, i) in used:
                    pass
                else:
                    used.append((i, second_dimension.m))
                    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_int(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_int(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)

@_string_unit_input
@ureg.wraps(("btuh", "ton_of_refrigeration"), ("person", "sqft", None), strict=False)
def cooling_load(ppl: int, sqft: float, external: bool = False):
    P_z = Q(ppl, "person")
    A_z = Q(sqft, "sqft")
    EL = Q("25 btuh / sqft") if external else Q("5 btuh / sqft")
    load = Q("450 btuh / person") * P_z + Q("1 W / sqft") * A_z + EL * A_z
    return load, load