Skip to content

Instantly share code, notes, and snippets.

@OrangeChannel

OrangeChannel/duct_sizer.py

Created December 26, 2022 05:49
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save OrangeChannel/f9aec0f0e0bd6371d20a68f7afe2479f to your computer and use it in GitHub Desktop.
Save OrangeChannel/f9aec0f0e0bd6371d20a68f7afe2479f to your computer and use it in GitHub Desktop.
Download ZIP
v8.1
Raw
duct_sizer.py
# 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, TypeAlias
from enum import Enum, auto
from matplotlib import pyplot as plt
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
import re
# endregion
# region --- Initializing Imports --------------------------------------------------------------------------------------
pd.options.display.max_columns = 20
F = TypeVar("F", bound=Callable)
ureg = UnitRegistry() # initialize a registry to hold the unit conversions
ureg.setup_matplotlib(True)
Q = ureg.Quantity # alias `Q` to the quantity constructor
QUANTITYLIKE: TypeAlias = Q | str | Real
width = Q(8.5, "inch")
margins = Q(1.75, "cm")
w = (width - 2 * margins) * 0.85
FIGSIZE = (w.to("inch").m, w.to("inch").m * 4 / 7)
# 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 SimpleRiser:
def __init__(self, riser_dict: dict[int, int], size_dict: Optional[dict[int, str]] = None):
"""
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.raw_riser_dict = riser_dict
self.raw_size_dict = size_dict or dict()
self._calc_cfms()
self._init_sizes()
def run(self, dims: Optional[list[int]] = None):
return run(self.cfms.m, dims)
def iupdate(self, riser_dict: dict[int, int], size_dict: Optional[dict[int, str]] = None):
self.raw_riser_dict |= riser_dict
if size_dict is None:
self.raw_size_dict = None
else:
self.raw_size_dict |= size_dict
self._calc_cfms()
self._init_sizes()
def update(self, riser_dict: dict[int, int], size_dict: Optional[dict[int, str]] = None):
self.iupdate(riser_dict, size_dict)
return self
def __add__(self, other: "SimpleRiser"):
if not isinstance(other, SimpleRiser):
raise NotImplementedError
new_riser_dict = dict()
for floor_number in self.raw_riser_dict:
if floor_number in other.raw_riser_dict:
new_riser_dict |= {
floor_number: self.raw_riser_dict[floor_number]
+ other.raw_riser_dict[floor_number]
}
return SimpleRiser(new_riser_dict)
def __getitem__(self, item):
return Q(self.cfms_dict[item], "cfm")
def _calc_cfms(self):
self.cfms_dict = dict()
for i, (floor, cfm) in enumerate(self.raw_riser_dict.items()):
self.cfms_dict[floor] = sum(list(self.raw_riser_dict.values())[:i + 1])
self.cfms = Q(list(self.cfms_dict.values()), "cfm")
def _init_sizes(self):
if self.raw_size_dict is None:
self.duct_dims = None
else:
self.duct_dims = dict()
for floor, dim in self.raw_size_dict.items():
self.duct_dims[floor] = DuctDimension(dim)
def verify(self, max_vel: str | Q = "1200 fpm"):
verification = dict()
for floor in self.raw_size_dict.keys():
if find_second_dimension(self[floor], self.duct_dims[floor][0], max_vel) <= self.duct_dims[floor][1]:
verification[floor] = True
else:
verification[floor] = False
return verification
def get_all_possible_dims(self, max_vel: str | Q = "1200 fpm"):
for (floor, cfm) in self.cfms_dict.items():
print(f"Floor {floor}:")
get_all_possible_dims(cfm, max_vel)
print(f"----------------")
def __repr__(self):
return f"<SimpleRiser({self.raw_riser_dict, self.raw_size_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.
TBD, quantify the cfms, make the cfms respect the floor numbering (dict), allow for use of a size_dict like Simple Riser
"""
self.connections_dict = connections_dict
self.branches = list(connections_dict.values())
self.total_connections = sum(self.branches)
self.cfm_per_connection = cfm_per_connection
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)
def __getitem__(self, item):
return self.cfms[item]
# endregion
# region --- Custom Unit Definitions -----------------------------------------------------------------------------------
# fmt: off
ureg.define("@alias sq_ft = sf = sqft")
ureg.define("@alias cu_ft = cf")
ureg.define("@alias Btu = btu")
ureg.define("cfm = cu_ft / min")
ureg.define("cfs = cu_ft / second")
ureg.define("fpm = ft / min")
ureg.define("@alias inch_H2O_39F = in_WC = inch_WC")
ureg.define("inch_WC_per_100_ft = inch_H2O_39F per (100 * ft)")
ureg.define("@alias inch_WC_per_100_ft = in_WC_per_100_ft")
ureg.define("person = [people]")
ureg.define("[occupant_density] = [people] / [area]")
ureg.define("ppl_per_1000_sqft = person per (1000 * sqft)")
ureg.define("btuh = Btu / hr")
ureg.define("mbh = 1000 * btuh")
ureg.define("percent = 0.01")
ureg.define("noise_criteria = [NC] = NC")
ureg.define("mcf = 1000 * cf")
ureg.define("brake_horsepower = hp = bhp")
# fmt: on
# endregion ------------------------------------------------------------------------------------------------------------
# region --- Add String Support for Future TUI/GUI ---------------------------------------------------------------------
def accept_strings(func: F) -> F:
@wraps(func)
def inner(*args: Any):
return func(*[(Q(x) if isinstance(x, str) else x) for x in args])
return inner
# endregion
# region --- Enums -----------------------------------------------------------------------------------------------------
class LinearDiffuserMounting(Enum):
SIDEWALL = auto()
CEILING = auto()
SIDEWALL = LinearDiffuserMounting.SIDEWALL
CEILING = LinearDiffuserMounting.CEILING
class PlenumBoxConnection(Enum):
CENTER = auto()
END = auto()
CENTER = PlenumBoxConnection.CENTER
END = PlenumBoxConnection.END
class AirFlowDirection(Enum):
SUPPLY = auto()
RETURN = auto()
SUPPLY = AirFlowDirection.SUPPLY
RETURN = AirFlowDirection.RETURN
# endregion
# region --- Architectual Dimension ------------------------------------------------------------------------------------
class ArchitectualLength:
def __init__(self, value1: int, value2: Optional[int] = None, value3: Optional[Real] = None, subprecision: Fraction = Fraction(1, 32)):
"""
Python representation of an architectural length quantity. Helpful for pretty-printing and quick inputting from
a REPL. Supports
Examples::
>>> q = Q(10.8188125, 'meter')
>>> s = r"16'-3 1/4\\"" # fix for dumb docstring bug in pycharm, string needs to be input correctly
>>> x, y, z = ArchitectualLength(3, 0, 0), ArchitectualLength(s), ArchitectualLength(q)
>>> x, y, z
(<ArchitectualLength(3, 0, 0)>, <ArchitectualLength(16, 3, 1/4)>, <ArchitectualLength(35, 5, 8444249301319921/9007199254740992)>)
>>> print(x, "\t\t", y, "\t\t", z)
3'-0 1/4" 16'-3 1/4" 35'-5 15/16"
>>> x.Q()
<Quantity(3.02083333, 'foot')>
"""
self.subprecision = subprecision
if (value2 is None) and (value3 is None):
if isinstance(value1, Q):
self._from_quantity(quant=value1)
elif isinstance(value1, str):
self._from_string(string=value1)
else:
self.feet = value1
self.inch = value2
self.subinch = Fraction(value3)
def as_quantity(self):
return Q(self.feet, "feet") + Q(self.inch, "inch") + Q(float(self.subinch), "inch")
Q = as_quantity
def __str__(self):
if (self.feet == 0) and (self.subinch == 0):
return f"{self.inch}\"" # 0" / 3"
elif (self.feet == 0) and (self.inch == 0):
return f"{self.subinch.limit_denominator(self.subprecision.denominator)}\"" # 1/4"
elif self.feet == 0:
return f"{self.inch} {self.subinch.limit_denominator(self.subprecision.denominator)}\""
elif (self.inch == 0) and (self.subinch == 0):
return f"{self.feet}'"
elif self.subinch == 0:
return f"{self.feet}'-{self.inch}\""
else:
return f"{self.feet}'-{self.inch} {self.subinch.limit_denominator(self.subprecision.denominator)}\""
def __repr__(self):
return f"<ArchitectualLength({self.feet}, {self.inch}, {self.subinch})>"
def _from_string(self, string: str):
p = re.compile(r"(?P<feet>\d+\')?-?(?P<inches>\d+[ \d/]*\")?")
m = p.search(string)
if m.group('feet'):
feet = int(m.group('feet')[:-1])
else:
feet = 0
if m.group('inches'):
_inches = m.group('inches').split(" ")
inches = int(_inches[0].replace('"', ''))
if len(_inches) == 2:
partial_inches = Fraction(_inches[1].replace('"', ''))
else:
partial_inches = 0
else:
inches = 0
partial_inches = 0
self.feet = feet
self.inch = inches
self.subinch = partial_inches
def _from_quantity(self, quant: Q):
feet = Q(math.floor(quant.to("ft").m), "ft")
self.feet = feet.m
quant -= feet
inches = Q(math.floor(quant.to("inch").m), "inch")
self.inch = inches.m
quant -= inches
self.subinch = Fraction(quant.to("inch").m)
# 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")
max_AR: Final[float] = 3.5
# endregion
# region --- Helper Functions ------------------------------------------------------------------------------------------
def _round_up_to_nearest_int(x: Real | Q, multiple: Real, abs_tol: float = 0.1) -> int:
"""
Rounds a numerical value, ``x``, up to the nearest ``multiple`` (integer). Will round down if ``x`` is only
``abs_tol`` greater than the closest multiple of ``multiple``.
:param x: Numerical value to be rounded up.
:param multiple: Integer multiple to round up to.
:param abs_tol: Absolute difference acceptable between ``x`` and the nearest multiple of ``multiple`` such that
``x`` is greater than this multiple and will be rounded down.
"""
if isinstance(x, Q):
if x.m < multiple:
return Q(multiple, x.u)
elif (x.m % multiple) < abs_tol:
return np.floor(x)
else:
return np.ceil(x / multiple) * multiple
else:
if x < multiple:
return multiple
elif (x % multiple) < abs_tol:
return math.floor(x)
else:
return math.ceil(x / multiple) * multiple
class DuctDimension:
def __init__(self, string: str):
string = string.lower()
L, W = string.split("x")
self.L = Q(f"{L} inch")
self.W = Q(f"{W} inch")
self.D_H = equiv_diameter_rect_duct(self.L, self.W)
def __repr__(self):
if self.L >= self.W:
return f"{self.L.m}x{self.W.m}"
else:
return f"{self.W.m}x{self.L.m}"
__str__ = __repr__
def __getitem__(self, item):
"""Makes it so you can use in functions that require both dimensions of the duct.
Example:
x = DuctDimension('8x6')
f(*x[:]) to call f with 8 inches and 6 inches as the first two arguments.
"""
return [self.L, self.W][item]
def stats(self, airflow: str | Q):
get_airflow_stats(self.L, self.W, airflow)
# endregion
# region --- Standalone Math Functions ---------------------------------------------------------------------------------
@accept_strings
@ureg.wraps("", ("m", "", "m"))
def Niazkar_friction_factor(
hydraulic_diameter: Q,
Reynolds_number: Q,
roughness_factor: Q = absolute_roughness,
) -> Q: # [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] (default set to ``absolute_roughness``)
:return: Niazkar friction factor [-]
"""
if math.isclose(hydraulic_diameter, 0, abs_tol=1e-9):
return 0
D = hydraulic_diameter
Re = Reynolds_number
ε = roughness_factor
p_1 = 4.5547
p_2 = 0.8784
A = -2 * math.log10(ε / (3.7 * D) + p_1 / pow(Re, p_2))
B = -2 * math.log10(ε / (3.7 * D) + 2.51 * A / Re)
C = -2 * math.log10(ε / (3.7 * D) + 2.51 * B / Re)
return pow(A - pow(B - A, 2) / (C - 2 * B + A), -2)
@accept_strings
@ureg.wraps("m", ("m^3 / s", "", "Pa / m", "kg / m^3"))
def hydraulic_diameter(
volumetric_flow_rate: Q,
friction_factor: Q,
head_loss: Q = max_head_loss,
fluid_density: Q = 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] (default set to ``max_head_loss``)
:param fluid_density: [kg/m^3] (default set to ``air_density``)
:return: Minimum pipe diameter [m]
"""
if math.isclose(volumetric_flow_rate, 0, abs_tol=1e-9):
return 0
f = friction_factor
ρ = fluid_density
H_L = head_loss
Q_dot = volumetric_flow_rate
return pow(
8 * pow(Q_dot, 2) * f * ρ / pow(np.pi, 2) / H_L,
1 / 5,
)
@accept_strings
@ureg.wraps("", ("m", "m/s", "m^2 / s"))
def Reynolds_number(
hydraulic_diameter: Q,
velocity: Q,
kinematic_viscosity: Q = 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 math.isclose(hydraulic_diameter, 0, abs_tol=1e-9):
return 0
D_h = hydraulic_diameter
V = velocity
ν = kinematic_viscosity
return D_h * V / ν
@accept_strings
@ureg.wraps("m/s", ("m^3 / s", "m^2"))
def fluid_velocity(volumetric_flow_rate: Q, orifice_area: Q) -> 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]
"""
if math.isclose(volumetric_flow_rate, 0, abs_tol=1e-9):
return 0
if math.isclose(orifice_area, 0, abs_tol=1e-9):
return 0
Q_dot = volumetric_flow_rate
A = orifice_area
return Q_dot / A
@accept_strings
@ureg.wraps("m^2", "m")
def pipe_area(internal_diameter: Q) -> Q:
"""
Calculate pipe area based on pipe diameter.
:param internal_diameter: [m]
:return: Pipe area [m^2]
"""
if math.isclose(internal_diameter, 0, abs_tol=1e-9):
return 0
D = internal_diameter
return math.pi / 4 * pow(D, 2)
@accept_strings
@ureg.wraps("inch", ("inch", "inch"))
def equiv_diameter_rect_duct(a: Q, b: Q) -> Q:
"""
Finds equivalent round duct given length and width of a rectangular duct (all in inches).
:param a: [inch]
:param b: [inch]
:return: [inch]
"""
if math.isclose(a, 0, abs_tol=1e-9):
return 0
if math.isclose(b, 0, abs_tol=1e-9):
return 0
return 1.30 * pow(a * b, 0.625) / pow(a + b, 0.250)
# endregion
# region --- Convenience Backend Wrappers ------------------------------------------------------------------------------
@accept_strings
def fluid_velocity_thru_pipe(volumetric_flow_rate: Q, internal_diameter: Q) -> Q:
"""
Wrapper for :py:func:fluid_velocity().
Calculate fluid velocity through a pipe given volumetric flow rate and pipe diameter.
"""
Q_dot = volumetric_flow_rate
ID = internal_diameter
A = pipe_area(ID)
return fluid_velocity(Q_dot, A)
@accept_strings
def fluid_velocity_thru_rect_duct_dims(volumetric_flow_rate: Q, a: Q, b: Q) -> Q:
"""
Wrapper for :py:func:fluid_velocity().
Calculate fluid velocity through a rectangular cross section.
"""
Q_dot = volumetric_flow_rate
D_h = equiv_diameter_rect_duct(a, b)
A = pipe_area(D_h)
return fluid_velocity(Q_dot, A)
@accept_strings
def find_hydraulic_diameter(
volumetric_flow_rate: Q,
head_loss: Q = max_head_loss,
fluid_density: Q = air_density,
roughness_factor: Q = absolute_roughness,
fluid_kinematic_viscosity: Q = air_kin_viscosity,
min_bound: Q = Q(1e-4, "m"),
max_bound: Q = Q(5, "m"),
) -> 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]
:param min_bound: [m]
:param max_bound: [m]
:return: Hydraulic diameter [m]
"""
if math.isclose(volumetric_flow_rate.m, 0, abs_tol=1e-9):
return 0
Q_dot = volumetric_flow_rate
H_L = head_loss
ρ = fluid_density
ε = roughness_factor
ν = fluid_kinematic_viscosity
min_bound = min_bound.to("m")
max_bound = max_bound.to("m")
@ureg.wraps(None, ("m^3/s", "kg/m^3", "Pa / m"))
def rhs_eq(volumetric_flow_rate, density, head_loss):
return 8 * pow(volumetric_flow_rate, 2) * density / pow(np.pi, 2) / head_loss
@ureg.wraps(None, ("m", ""))
def lhs_eq(hydraulic_diameter, friction_factor):
return pow(hydraulic_diameter, 5) / friction_factor
rhs = rhs_eq(Q_dot, ρ, H_L)
def f(x):
"""Make sure `x` is in meters"""
x = Q(x, "m")
A = pipe_area(x)
V = fluid_velocity(Q_dot, A)
Re = Reynolds_number(x, V, ν)
N_f = Niazkar_friction_factor(x, Re, ε)
lhs = lhs_eq(x, N_f)
return abs(rhs - lhs)
return Q(minimize_scalar(f, method="bounded", bounds=(min_bound.m, max_bound.m)).x, "m")
@accept_strings
def find_second_dimension(
volumetric_flow_rate: Q,
first_dimension: Q,
max_vel: Q = max_vel,
head_loss: Q = max_head_loss,
fluid_density: Q = air_density,
roughness_factor: Q = absolute_roughness,
fluid_kinematic_viscosity: Q = air_kin_viscosity,
min_bound: Q = Q(1, "inch"),
max_bound: Q = Q(200, "inch"),
) -> 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: [m/s]
:param head_loss: [Pa/m]
:param fluid_density: [kg/m^3]
:param roughness_factor: [m]
:param fluid_kinematic_viscosity: [m^2/s]
:param: min_bound: [inch]
:param: max_bound: [inch]
:return: Other dimension in a rectangular duct [inch]
"""
if math.isclose(volumetric_flow_rate.m, 0, abs_tol=1e-9):
return 0
if math.isclose(first_dimension.m, 0, abs_tol=1e-9):
return 0
Q_dot = volumetric_flow_rate
a = first_dimension
V_max = max_vel
H_L = head_loss
ρ = fluid_density
ε = roughness_factor
ν = fluid_kinematic_viscosity
min_bound = min_bound.to("inch")
max_bound = max_bound.to("inch")
D_h = find_hydraulic_diameter(Q_dot, H_L, ρ, ε, ν).to("inch")
def f1(x): # minimizes for hydraulic dia
"""Make sure `x` is in inches"""
x = Q(x, "inch")
return abs(D_h.m - equiv_diameter_rect_duct(a, x).m)
def f2(x): # minimizes for velocity
"""Make sure `x` is inches"""
x = Q(x, "inch")
V = fluid_velocity_thru_rect_duct_dims(Q_dot, a, x)
return abs(V.m - max_vel.to("m/s").m)
second_dimension = Q(minimize_scalar(f1, method="bounded", bounds=(min_bound.m, max_bound.m)).x, "inch")
f1_air_vel = fluid_velocity_thru_rect_duct_dims(Q_dot, a, second_dimension)
if f1_air_vel <= V_max:
return second_dimension
else:
return Q(minimize_scalar(f2, method="bounded", bounds=(min_bound.m, max_bound.m)).x, "inch")
# endregion
# region --- User Functions --------------------------------------------------------------------------------------------
@accept_strings
@ureg.wraps(None, ("cfm", "fpm", "in_WC_per_100_ft", "inch", ""), strict=False)
def get_all_possible_dims(
volumetric_flow_rate: QUANTITYLIKE,
max_vel: QUANTITYLIKE = max_vel,
head_loss: QUANTITYLIKE = max_head_loss,
max_dim: QUANTITYLIKE = Q("48 inch"),
max_aspect_ratio: Real = max_AR,
):
if math.isclose(volumetric_flow_rate, 0, abs_tol=1e-9):
print("0 x 0")
return
Q_dot = Q(volumetric_flow_rate, "cfm")
V_max = Q(max_vel, "fpm")
H_L = Q(head_loss, "in_WC_per_100_ft")
used = []
for i in ([5] + list(range(6, max_dim + 1))[::2]) * ureg.inch:
second_dimension = find_second_dimension(Q_dot, i, V_max, H_L)
second_dimension = _round_up_to_nearest_int(second_dimension, 2)
if second_dimension >= i:
AR = (second_dimension / i).to("")
else:
AR = (i / second_dimension).to("")
if AR <= max_aspect_ratio:
if (int(second_dimension.m), int(i.m)) in used:
pass
elif int(second_dimension.m) in [u[0] for u in used]:
pass
else:
used.append((int(i.m), int(second_dimension.m)))
if AR == 1:
print(int(i.m), "x", int(second_dimension.m), " *")
else:
print(int(i.m), "x", int(second_dimension.m))
@accept_strings
@ureg.wraps(None, ("cfm", "in_WC_per_100_ft", None), strict=False)
def get_round_duct_size(
volumetric_flow_rate: QUANTITYLIKE,
head_loss: QUANTITYLIKE = max_head_loss,
print_: bool = True,
):
Q_dot = Q(volumetric_flow_rate, "cfm")
H_L = Q(head_loss, "in_WC_per_100_ft")
D_h = find_hydraulic_diameter(Q_dot, H_L).to("inch")
D_h = _round_up_to_nearest_int(D_h, 1)
if print_:
print(D_h)
else:
return D_h
@accept_strings
@ureg.wraps(None, ("cfm", "inch", "fpm", "in_WC_per_100_ft", None), strict=False)
def get_second_dimension(
volumetric_flow_rate: QUANTITYLIKE,
first_dimension: QUANTITYLIKE,
max_vel: QUANTITYLIKE = max_vel,
head_loss: QUANTITYLIKE = max_head_loss,
print_: bool = True,
):
Q_dot = Q(volumetric_flow_rate, "cfm")
a = Q(first_dimension, "inch")
V_max = Q(max_vel, "fpm")
H_L = Q(head_loss, "in_WC_per_100_ft")
b = find_second_dimension(Q_dot, a, V_max, H_L)
b = _round_up_to_nearest_int(b, 2)
if print_:
print(b)
else:
return b
@ureg.wraps(
None,
(
"sqft",
None,
"person",
"sqft / person",
None,
None,
None,
None,
"cfm"
),
strict=False,
)
def minimum_ventilation_rates(
zone_floor_area: Q,
occupancy_classification: str = 'Multipurpose assembly',
exact_occupancy: Q = Q("0 person"),
exact_area_per_person: Q = Q("0 sqft / person"),
zone_air_distribution_effectiveness: float = 0.8,
bedrooms: int = 0,
use_max_of_default: bool = False,
code_year: int = 2014,
minimum_flow_per_space: Q = Q("20 cfm")
):
# re-assign units
A_z = Q(zone_floor_area, "sqft")
P_z = Q(exact_occupancy, "person")
exact_density = Q(exact_area_per_person, "sqft / person")
min_flow = Q(minimum_flow_per_space, "cfm")
# import the 2014 data
table_403 = pd.read_csv(r'C:\Users\Dave\PycharmProjects\Ventrop\table_403_3_minimum_ventilation_rates.csv')
table_403 = table_403.set_index('Occupancy Classification')
# get the row based on room classifications
classification_row = table_403.loc[[occupancy_classification]]
R_p = Q(float(classification_row['R_p [cfm/person]']), "cfm/person")
R_a = Q(float(classification_row['R_a [cfm/sqft]']), "cfm/sqft")
default_occupant_density = Q(float(classification_row['Occupant Density [ppl_per_1000_sqft]']), "ppl_per_1000_sqft")
exhaust_airflow_rate = Q(float(classification_row['Exhaust Airflow [cfm/sqft]']), "cfm/sqft")
# occupancy calcs
if exact_density > Q(0, "sqft / person"):
P_z = np.ceil((A_z / exact_density).to("person"))
if use_max_of_default:
if occupancy_classification == 'Living areas':
P_z = max(P_z, Q(bedrooms + 1, "person"))
else:
P_z = max(P_z, np.ceil((default_occupant_density * A_z).to("person")))
V_bz = _round_up_to_nearest_int((R_p * P_z + R_a * A_z).to("cfm") / zone_air_distribution_effectiveness, 5)
V_bz = max(min_flow, V_bz)
exhaust_airflow = (exhaust_airflow_rate * A_z).to("cfm")
if exhaust_airflow > Q(0, "cfm"):
exhaust_airflow = max(min_flow, exhaust_airflow)
print(A_z, occupancy_classification)
if exhaust_airflow == Q(0, "cfm"):
print(f"{V_bz} OA/RA, {P_z}")
else:
exhaust_airflow = _round_up_to_nearest_int(exhaust_airflow, 5)
print(f"{V_bz} OA, {exhaust_airflow} EX, {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)
return df
@accept_strings
@ureg.wraps(("btuh", "ton_of_refrigeration"), ("person", "sqft", None), strict=False)
def cooling_load_estimation(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
@ureg.wraps("NC", ("feet", None), strict=False)
def _NC_correction(
plenum_box_length: QUANTITYLIKE, airflow_direction: AirFlowDirection = SUPPLY
) -> Q:
L = Q(plenum_box_length, "ft")
df_NC_correction = pd.DataFrame()
df_NC_correction["Length [ft]"] = [2, 4, 6, 8, 10]
df_NC_correction["Supply [NC]"] = [-2, 0, 2, 3, 5]
df_NC_correction["Return [NC]"] = [0, 3, 5, 6, 8]
match airflow_direction:
case AirFlowDirection.SUPPLY:
return np.interp(
L,
np.asarray(df_NC_correction["Length [ft]"]) * ureg.ft,
np.asarray(df_NC_correction["Supply [NC]"]) * ureg.NC,
)
case AirFlowDirection.RETURN:
return np.interp(
L,
np.asarray(df_NC_correction["Length [ft]"]) * ureg.ft,
np.asarray(df_NC_correction["Return [NC]"]) * ureg.NC,
)
@ureg.wraps(None, "feet", strict=False)
def _throw_correction_multiplier(plenum_box_length: QUANTITYLIKE) -> float:
L = Q(plenum_box_length, "ft")
df_throw_correction = pd.DataFrame()
df_throw_correction["Length [ft]"] = [2, 4, 8, 10, 12]
df_throw_correction["Throw Correction"] = [0.72, 1, 1.5, 1.7, 1.8]
return np.interp(
L,
np.asarray(df_throw_correction["Length [ft]"]) * ureg.ft,
df_throw_correction["Throw Correction"],
)
@accept_strings
@ureg.wraps(None, ("cfm", "inch", None, None, "NC", None), strict=False)
def linear_air_sizer(
airflow: QUANTITYLIKE,
incoming_duct_height: QUANTITYLIKE,
connection_type: PlenumBoxConnection,
mounting_type: LinearDiffuserMounting,
max_acceptable_noise: QUANTITYLIKE = Q(18, "NC"),
slots: int = 1,
):
Q_dot = Q(airflow, "cfm")
DH = Q(incoming_duct_height, "inch")
match mounting_type:
case LinearDiffuserMounting.SIDEWALL:
match slots:
case 1:
plenum_box_height = max(DH + Q("2 inch"), Q("5.75 inch"))
case 2:
plenum_box_height = max(
DH + Q("2 inch"), Q("11.25 inch")
) # see FBP submittal
match connection_type:
case PlenumBoxConnection.CENTER:
plenum_box_width = get_second_dimension(
Q_dot / 2,
plenum_box_height,
Q("1000 fpm"),
"0.1 in_WC_per_100_ft",
False,
)
case PlenumBoxConnection.END:
plenum_box_width = get_second_dimension(
Q_dot,
plenum_box_height,
Q("1000 fpm"),
"0.1 in_WC_per_100_ft",
False,
)
case LinearDiffuserMounting.CEILING:
plenum_box_height = DH + Q("2 inch")
match slots:
case 1:
min_width = Q("5.75 inch")
case 2:
min_width = Q("11.25 inch")
match connection_type:
case PlenumBoxConnection.CENTER:
plenum_box_width = max(
min_width,
get_second_dimension(
Q_dot / 2,
plenum_box_height,
Q("1000 fpm"),
"0.1 in_WC_per_100_ft",
False,
),
)
case PlenumBoxConnection.END:
plenum_box_width = max(
min_width,
get_second_dimension(
Q_dot,
plenum_box_height,
Q("1000 fpm"),
"0.1 in_WC_per_100_ft",
False,
),
)
max_NC = Q(max_acceptable_noise, "NC")
FL_25_specs = dict()
FL_25_specs[CEILING] = dict()
FL_25_specs[SIDEWALL] = dict()
FL_25_specs[CEILING][1] = dict()
FL_25_specs[CEILING][2] = dict()
FL_25_specs[SIDEWALL][1] = dict()
FL_25_specs[SIDEWALL][2] = dict()
FL_25_specs[CEILING][1]["Airflow [cfm/ft]"] = Q([30, 55, 80, 105, 130, 155, 180], "cfm/ft")
FL_25_specs[CEILING][1]["Noise Criteria [NC]"] = Q([10, 11, 23, 31, 38, 43, 48], "NC")
FL_25_specs[CEILING][1]["Throw (L) [ft]"] = Q([4, 10, 14, 17, 19, 21, 22], "ft")
FL_25_specs[CEILING][1]["Throw (M) [ft]"] = Q([8, 15, 18, 21, 23, 25, 27], "ft")
FL_25_specs[CEILING][1]["Throw (H) [ft]"] = Q([16, 21, 26, 29, 33, 36, 39], "ft")
FL_25_specs[CEILING][2]["Airflow [cfm/ft]"] = Q([60, 105, 150, 195, 240, 285, 330], "cfm/ft")
FL_25_specs[CEILING][2]["Noise Criteria [NC]"] = Q([10, 13, 24, 32, 38, 44, 48], "NC")
FL_25_specs[CEILING][2]["Throw (L) [ft]"] = Q([6, 13, 19, 23, 26, 28, 30], "ft")
FL_25_specs[CEILING][2]["Throw (M) [ft]"] = Q([11, 20, 25, 28, 32, 34, 37], "ft")
FL_25_specs[CEILING][2]["Throw (H) [ft]"] = Q([22, 29, 35, 40, 45, 49, 52], "ft")
FL_25_specs[SIDEWALL][1]["Airflow [cfm/ft]"] = Q([95, 150, 205, 260, 315, 370], "cfm/ft")
FL_25_specs[SIDEWALL][1]["Noise Criteria [NC]"] = Q([10, 18, 25, 30, 35, 39], "NC")
FL_25_specs[SIDEWALL][1]["Throw (L) [ft]"] = Q([6, 9, 12, 15, 18, 22], "ft")
FL_25_specs[SIDEWALL][1]["Throw (M) [ft]"] = Q([8, 13, 18, 23, 28, 33], "ft")
FL_25_specs[SIDEWALL][1]["Throw (H) [ft]"] = Q([17, 26, 35, 39, 43, 46], "ft")
FL_25_specs[SIDEWALL][2]["Airflow [cfm/ft]"] = Q([80, 190, 300, 410, 520, 630, 740], "cfm/ft")
FL_25_specs[SIDEWALL][2]["Noise Criteria [NC]"] = Q([10, 11, 22, 28, 34, 38, 42], "NC")
FL_25_specs[SIDEWALL][2]["Throw (L) [ft]"] = Q([2, 8, 12, 17, 22, 26, 31], "ft")
FL_25_specs[SIDEWALL][2]["Throw (M) [ft]"] = Q([4, 12, 19, 26, 32, 39, 46], "ft")
FL_25_specs[SIDEWALL][2]["Throw (H) [ft]"] = Q([10, 24, 37, 49, 55, 60, 66], "ft")
flow_per_ft_adj = lambda x: np.interp(max_NC - x,
FL_25_specs[mounting_type][slots]["Noise Criteria [NC]"],
FL_25_specs[mounting_type][slots]["Airflow [cfm/ft]"],
)
max_airflow_per_length = flow_per_ft_adj(Q("0 NC"))
prev_length = Q_dot / max_airflow_per_length
NC_adjustment = _NC_correction(prev_length)
new_max_airflow_per_length = flow_per_ft_adj(NC_adjustment)
new_length = Q_dot / new_max_airflow_per_length
while abs(new_length - prev_length) > Q("1/2 inch"):
prev_length = new_length
NC_adjustment = _NC_correction(prev_length)
new_max_airflow_per_length = flow_per_ft_adj(NC_adjustment)
new_length = Q_dot / new_max_airflow_per_length
plenum_box_length = Q(_round_up_to_nearest_int(new_length.to("inch").m, 4), "inch")
print(
f"Plenum box dimensions (LxWxH): {plenum_box_length} x {plenum_box_width} x {plenum_box_height}"
)
LD_supply_length = Q(
_round_up_to_nearest_int(plenum_box_length.m + 2, 6), "inch"
).to("ft")
max_airflow_per_length = flow_per_ft_adj(Q("0 NC"))
prev_length = Q_dot / max_airflow_per_length
NC_adjustment = _NC_correction(prev_length, airflow_direction=RETURN)
new_max_airflow_per_length = flow_per_ft_adj(NC_adjustment)
new_length = Q_dot / new_max_airflow_per_length
while abs(new_length - prev_length) > Q("1/2 inch"):
prev_length = new_length
NC_adjustment = _NC_correction(prev_length, airflow_direction=RETURN)
new_max_airflow_per_length = flow_per_ft_adj(NC_adjustment)
new_length = Q_dot / new_max_airflow_per_length
LD_return_length = Q(
_round_up_to_nearest_int(new_length.to("inch").m, 6), "inch"
).to("ft")
total_LD_length = LD_supply_length + LD_return_length
print(
f"Linear diffuser length: ({LD_supply_length} SUPPLY + {LD_return_length} RETURN) = {total_LD_length}"
)
print(f"If combined LD: LD-{slots} // {ArchitectualLength(total_LD_length)}")
print(
f"If separate LD: Supply = LD-{slots} // {ArchitectualLength(LD_supply_length)}\tReturn = LD-{slots} // {ArchitectualLength(LD_return_length)}"
)
throw_dist = dict()
for k in "LMH":
throw_dist[k] = np.interp(
Q_dot / plenum_box_length.to("ft"),
FL_25_specs[mounting_type][slots]["Airflow [cfm/ft]"],
FL_25_specs[mounting_type][slots][f"Throw ({k}) [ft]"] * _throw_correction_multiplier(plenum_box_length=plenum_box_length),
)
print(
f"Throw distance (L, M, H): {ArchitectualLength(throw_dist['L'])}, {ArchitectualLength(throw_dist['M'])}, {ArchitectualLength(throw_dist['H'])}"
)
@accept_strings
@ureg.wraps(None, ("cfm", "inch", "", "fpm"), strict=False)
def louver_sizing(airflow: Q, first_dimension: Q, free_area: Q = Q("50 percent"), airspeed: Q = Q("600 fpm")):
Q_dot = Q(airflow, "cfm")
a = Q(first_dimension, "inch")
v = Q(airspeed, "fpm")
f = free_area
b = (Q_dot / f / a / v).to("inch")
return np.ceil(b), np.round((a*b).to("sqft"), 1)
@accept_strings
@ureg.wraps(None, ("sqft", "ft", "lbs"), strict=False)
def refrigerant_leak_check_R410A(floor_area: Q, floor_to_ceiling: Q, total_charge: Q):
"""
>>> refrigerant_leak_check_R410A("315 sqft", ArchitectualLength("10'-2\\"").Q(), "76 lbs")
# SAFE
"""
max_conc = Q("26 lb / mcf")
room_vol = (floor_area * floor_to_ceiling).to("mcf")
charge_per_vol = (total_charge / room_vol).to("lb / mcf")
if charge_per_vol <= max_conc:
print("SAFE")
else:
print("NOT SAFE")
@accept_strings
@ureg.wraps(None, ("inch", "inch", "cfm"), strict=False)
def get_airflow_stats(a: Q, b: Q, airflow: Q):
a = Q(a, "inch")
b = Q(b, "inch")
airflow = Q(airflow, "cfm")
D_e = equiv_diameter_rect_duct(a, b).to("inch")
flow_area = (a * b).to("ft^2")
fluid_vel = fluid_velocity_thru_rect_duct_dims(airflow, a, b).to("fpm")
Re = Reynolds_number(D_e, fluid_vel)
f = Niazkar_friction_factor(D_e, Re)
vel_pressure = (air_density * pow(fluid_vel, 2) / 2).to("in_WC")
head_loss = (pow(fluid_vel, 2) * f * air_density / 2 / D_e).to("in_WC_per_100_ft")
print(f"For {airflow:.0f~P} thru a {a} x {b} duct:"
f"\n\tEquivalent diameter = {D_e:.3f~P}"
f"\n\tFlow area = {flow_area:.4f~P}"
f"\n\tFluid velocity = {fluid_vel:,.0f~P}"
f"\n\tReynold's number = {Re:,.0f~P}"
f"\n\tFriction factor = {f:.5f~P}"
f"\n\tVelocity pressure = {vel_pressure:.4f~P}"
f"\n\tHead loss = {head_loss:.4f~P}")
@ZeZeZe99
Copy link

ZeZeZe99 commented Apr 4, 2024

Hi Dave, I'm wondering what is the equation you use for calculating the pipe diameter in line 450-453. I cannot find a formula that is similar to this one. Thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment