"""Chemical Engineering Design Library (ChEDL). Utilities for process modeling. Copyright (C) 2016, 2017, 2018, 2019, 2020 Caleb Bell Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. This module contains equations for modeling control valves subject to gas or liquid flow. For reporting bugs, adding feature requests, or submitting pull requests, please use the `GitHub issue tracker `_ or contact the author at Caleb.Andrew.Bell@gmail.com. .. contents:: :local: Sizing Functions ---------------- .. autofunction:: size_control_valve_l .. autofunction:: size_control_valve_g Intermediary Sizing Calculations -------------------------------- .. autofunction:: FF_critical_pressure_ratio_l .. autofunction:: is_choked_turbulent_l .. autofunction:: is_choked_turbulent_g .. autofunction:: Reynolds_valve .. autofunction:: Reynolds_factor .. autofunction:: loss_coefficient_piping .. autofunction:: control_valve_choke_P_l .. autofunction:: control_valve_choke_P_g .. autofunction:: convert_flow_coefficient .. autofunction:: cavitation_index Representative Control Valve Curves ----------------------------------- .. autofunction:: Cv_char_linear .. autofunction:: Cv_char_quick_opening .. autofunction:: Cv_char_equal_percentage Noise Generated by Control Valves --------------------------------- .. autofunction:: control_valve_noise_l_2015 .. autofunction:: control_valve_noise_g_2011 """ from math import exp, log, log10, pi, sqrt from fluids.constants import R, gallon, ln_10, ln_10_inv, minute, psi from fluids.fittings import Cv_to_Kv, Kv_to_Cv from fluids.numerics import implementation_optimize_tck, interp, splev __all__ = ['size_control_valve_l', 'size_control_valve_g', 'cavitation_index', 'FF_critical_pressure_ratio_l', 'is_choked_turbulent_l', 'is_choked_turbulent_g', 'Reynolds_valve', 'loss_coefficient_piping', 'Reynolds_factor', 'Cv_char_quick_opening', 'Cv_char_linear', 'Cv_char_equal_percentage', 'convert_flow_coefficient', 'control_valve_choke_P_l', 'control_valve_choke_P_g', 'control_valve_noise_l_2015', 'control_valve_noise_g_2011'] N1 = 0.1 # m^3/hr, kPa N2 = 1.6E-3 # mm N4 = 7.07E-2 # m^3/hr, m^2/s N5 = 1.8E-3 # mm N6 = 3.16 # kg/hr, kPa, kg/m^3 N7 = 4.82 # m^3/hr kPa K N8 = 1.10 # kPa kg/hr K #N9 = 2.60E1 # m^3/hr kPa K at 15 deg C N9 = 2.46E1 # m^3/hr kPa K at 0 deg C N18 = 8.65E-1 # mm N19 = 2.5 # mm #N22 = 1.84E1 # m^3/hr kPa K at 15 deg C N27 = 7.75E-1 # kg/hr kPa K at 0 deg C N32 = 1.4E2 # mm rho0 = 999.10329075702327 # Water at 288.15 K def cavitation_index(P1, P2, Psat): r'''Calculates the cavitation index of a valve with upstream and downstream absolute pressures `P1` and `P2` for a fluid with a vapor pressure `Psat`. .. math:: \sigma = \frac{P_1 - P_{sat}}{P_1 - P_2} Parameters ---------- P1 : float Absolute pressure upstream of the valve [Pa] P2 : float Absolute pressure downstream of the valve [Pa] Psat : float Saturation pressure of the liquid at inlet temperature [Pa] Returns ------- sigma : float Cavitation index of the valve [-] Notes ----- Larger values are safer. Models for adjusting cavitation indexes provided by the manufacturer to the user's conditions are available, making use of scaling the pressure differences and size differences. Values can be calculated for incipient cavitation, constant cavitation, maximum vibration cavitation, incipient damage, and choking cavitation. Has also been defined as: .. math:: \sigma = \frac{P_2 - P_{sat}}{P_1 - P_2} Another definition and notation series is: .. math:: K = xF = \frac{1}{\sigma} = \frac{P_1 - P_2}{P_1 - P_{sat}} Examples -------- >>> cavitation_index(1E6, 8E5, 2E5) 4.0 References ---------- .. [1] ISA. "RP75.23 Considerations for Evaluating Control Valve Cavitation." 1995. ''' return (P1 - Psat)/(P1 - P2) def FF_critical_pressure_ratio_l(Psat, Pc): r'''Calculates FF, the liquid critical pressure ratio factor, for use in IEC 60534 liquid valve sizing calculations. .. math:: F_F = 0.96 - 0.28\sqrt{\frac{P_{sat}}{P_c}} Parameters ---------- Psat : float Saturation pressure of the liquid at inlet temperature [Pa] Pc : float Critical pressure of the liquid [Pa] Returns ------- FF : float Liquid critical pressure ratio factor [-] Examples -------- From [1]_, matching example. >>> FF_critical_pressure_ratio_l(70100.0, 22120000.0) 0.9442375225233299 References ---------- .. [1] IEC 60534-2-1 / ISA-75.01.01-2007 ''' return 0.96 - 0.28*sqrt(Psat/Pc) def control_valve_choke_P_l(Psat, Pc, FL, P1=None, P2=None, disp=True): r'''Calculates either the upstream or downstream pressure at which choked flow though a liquid control valve occurs, given either a set upstream or downstream pressure. Implements an analytical solution of the needed equations from the full function :py:func:`~.size_control_valve_l`. For some pressures, no choked flow is possible; for choked flow to occur the direction if flow must be reversed. If `disp` is True, an exception will be raised for these conditions. .. math:: P_1 = \frac{F_{F} F_{L}^{2} P_{sat} - P_{2}}{F_{L}^{2} - 1} .. math:: P_2 = F_{F} F_{L}^{2} P_{sat} - F_{L}^{2} P_{1} + P_{1} Parameters ---------- Psat : float Saturation pressure of the liquid at inlet temperature [Pa] Pc : float Critical pressure of the liquid [Pa] FL : float, optional Liquid pressure recovery factor of a control valve without attached fittings [-] P1 : float, optional Absolute pressure upstream of the valve [Pa] P2 : float, optional Absolute pressure downstream of the valve [Pa] disp : bool, optional Whether or not to raise an exception on flow reversal, [-] Returns ------- P_choke : float Pressure at which a choke occurs in the liquid valve [Pa] Notes ----- Extremely cheap to compute. Examples -------- >>> control_valve_choke_P_l(69682.89291024722, 22048320.0, 0.6, 680000.0) 458887.5306077305 >>> control_valve_choke_P_l(69682.89291024722, 22048320.0, 0.6, P2=458887.5306077305) 680000.0 ''' FF = 0.96 - 0.28*sqrt(Psat/Pc) #FF_critical_pressure_ratio_l(Psat=Psat, Pc=Pc) Pmin_absolute = FF*Psat if P2 is None: ans = P2 = FF*FL*FL*Psat - FL*FL*P1 + P1 elif P1 is None: ans = P1 = (FF*FL*FL*Psat - P2)/(FL*FL - 1.0) else: raise ValueError('Either P1 or P2 needs to be specified') if P2 > P1 and disp: raise ValueError('Specified P1 is too low for choking to occur ' 'at any downstream pressure; minimum ' 'upstream pressure for choking to be possible ' f'is {Pmin_absolute:g} Pa.') return ans def control_valve_choke_P_g(xT, gamma, P1=None, P2=None): r'''Calculates either the upstream or downstream pressure at which choked flow though a gas control valve occurs, given either a set upstream or downstream pressure. Implements an analytical solution of the needed equations from the full function :py:func:`~.size_control_valve_g`. A singularity arises as `xT` goes to 1 and `gamma` goes to 1.4. .. math:: P_1 = - \frac{7 P_{2}}{5 \gamma x_T - 7} .. math:: P_2 = \frac{P_{1}}{7} \left(- 5 \gamma x_T + 7\right) Parameters ---------- xT : float, optional Pressure difference ratio factor of a valve without fittings at choked flow [-] gamma : float Specific heat capacity ratio [-] P1 : float, optional Absolute pressure upstream of the valve [Pa] P2 : float, optional Absolute pressure downstream of the valve [Pa] Returns ------- P_choke : float Pressure at which a choke occurs in the gas valve [Pa] Notes ----- Extremely cheap to compute. Examples -------- >>> control_valve_choke_P_g(1.0, 1.3, 1E5) 7142.857142857143 >>> control_valve_choke_P_g(1.0, 1.3, P2=7142.857142857143) 100000.0 ''' if P2 is None: ans = P2 = P1*(-5.0*gamma*xT + 7.0)/7.0 elif P1 is None: ans = P1 = -7.0*P2/(5.0*gamma*xT - 7.0) else: raise ValueError('Either P1 or P2 needs to be specified') return ans def is_choked_turbulent_l(dP, P1, Psat, FF, FL=None, FLP=None, FP=None): r'''Calculates if a liquid flow in IEC 60534 calculations is critical or not, for use in IEC 60534 liquid valve sizing calculations. Either FL may be provided or FLP and FP, depending on the calculation process. .. math:: \Delta P > F_L^2(P_1 - F_F P_{sat}) .. math:: \Delta P >= \left(\frac{F_{LP}}{F_P}\right)^2(P_1 - F_F P_{sat}) Parameters ---------- dP : float Differential pressure across the valve, with reducer/expanders [Pa] P1 : float Pressure of the fluid before the valve and reducers/expanders [Pa] Psat : float Saturation pressure of the fluid at inlet temperature [Pa] FF : float Liquid critical pressure ratio factor [-] FL : float, optional Liquid pressure recovery factor of a control valve without attached fittings [-] FLP : float, optional Combined liquid pressure recovery factor with piping geometry factor, for a control valve with attached fittings [-] FP : float, optional Piping geometry factor [-] Returns ------- choked : bool Whether or not the flow is choked [-] Examples -------- >>> is_choked_turbulent_l(460.0, 680.0, 70.1, 0.94, 0.9) False >>> is_choked_turbulent_l(460.0, 680.0, 70.1, 0.94, 0.6) True References ---------- .. [1] IEC 60534-2-1 / ISA-75.01.01-2007 ''' if FLP and FP: return dP >= FLP*FLP/(FP*FP)*(P1-FF*Psat) elif FL: return dP >= FL*FL*(P1-FF*Psat) else: raise ValueError('Either (FLP and FP) or FL is needed') def is_choked_turbulent_g(x, Fgamma, xT=None, xTP=None): r'''Calculates if a gas flow in IEC 60534 calculations is critical or not, for use in IEC 60534 gas valve sizing calculations. Either xT or xTP must be provided, depending on the calculation process. .. math:: x \ge F_\gamma x_T .. math:: x \ge F_\gamma x_{TP} Parameters ---------- x : float Differential pressure over inlet pressure, [-] Fgamma : float Specific heat ratio factor [-] xT : float, optional Pressure difference ratio factor of a valve without fittings at choked flow [-] xTP : float Pressure difference ratio factor of a valve with fittings at choked flow [-] Returns ------- choked : bool Whether or not the flow is choked [-] Examples -------- Example 3, compressible flow, non-choked with attached fittings: >>> is_choked_turbulent_g(0.544, 0.929, 0.6) False >>> is_choked_turbulent_g(0.544, 0.929, xTP=0.625) False References ---------- .. [1] IEC 60534-2-1 / ISA-75.01.01-2007 ''' if xT: return x >= Fgamma*xT elif xTP: return x >= Fgamma*xTP else: raise ValueError('Either xT or xTP is needed') def Reynolds_valve(nu, Q, D1, FL, Fd, C): r'''Calculates Reynolds number of a control valve for a liquid or gas flowing through it at a specified Q, for a specified D1, FL, Fd, C, and with kinematic viscosity `nu` according to IEC 60534 calculations. .. math:: Re_v = \frac{N_4 F_d Q}{\nu \sqrt{C F_L}}\left(\frac{F_L^2 C^2} {N_2D^4} +1\right)^{1/4} Parameters ---------- nu : float Kinematic viscosity, [m^2/s] Q : float Volumetric flow rate of the fluid [m^3/s] D1 : float Diameter of the pipe before the valve [m] FL : float, optional Liquid pressure recovery factor of a control valve without attached fittings [] Fd : float Valve style modifier [-] C : float Metric Kv valve flow coefficient (flow rate of water at a pressure drop of 1 bar) [m^3/hr] Returns ------- Rev : float Valve reynolds number [-] Examples -------- >>> Reynolds_valve(3.26e-07, 360, 150.0, 0.9, 0.46, 165) 2966984.7525455453 References ---------- .. [1] IEC 60534-2-1 / ISA-75.01.01-2007 ''' return N4*Fd*Q/nu*1.0/sqrt(C*FL)*sqrt(sqrt(FL*FL*C*C/N2*D1**-4.0 + 1.0)) def loss_coefficient_piping(d, D1=None, D2=None): r'''Calculates the sum of loss coefficients from possible inlet/outlet reducers/expanders around a control valve according to IEC 60534 calculations. .. math:: \Sigma \xi = \xi_1 + \xi_2 + \xi_{B1} - \xi_{B2} .. math:: \xi_1 = 0.5\left[1 -\left(\frac{d}{D_1}\right)^2\right]^2 .. math:: \xi_2 = 1.0\left[1 -\left(\frac{d}{D_2}\right)^2\right]^2 .. math:: \xi_{B1} = 1 - \left(\frac{d}{D_1}\right)^4 .. math:: \xi_{B2} = 1 - \left(\frac{d}{D_2}\right)^4 Parameters ---------- d : float Diameter of the valve [m] D1 : float Diameter of the pipe before the valve [m] D2 : float Diameter of the pipe after the valve [m] Returns ------- loss : float Sum of the four loss coefficients [-] Examples -------- In example 3, non-choked compressible flow with fittings: >>> loss_coefficient_piping(0.05, 0.08, 0.1) 0.6580810546875 References ---------- .. [1] IEC 60534-2-1 / ISA-75.01.01-2007 ''' loss = 0. if D1: dr = d/D1 dr2 = dr*dr loss += 1. - dr2*dr2 # Inlet flow energy loss += 0.5*(1. - dr2)*(1.0 - dr2) # Inlet reducer if D2: dr = d/D2 dr2 = dr*dr loss += 1.0*(1. - dr2)*(1.0 - dr2) # Outlet reducer (expander) loss -= 1. - dr2*dr2 # Outlet flow energy return loss def Reynolds_factor(FL, C, d, Rev, full_trim=True): r'''Calculates the Reynolds number factor `FR` for a valve with a Reynolds number `Rev`, diameter `d`, flow coefficient `C`, liquid pressure recovery factor `FL`, and with either full or reduced trim, all according to IEC 60534 calculations. If full trim: .. math:: F_{R,1a} = 1 + \left(\frac{0.33F_L^{0.5}}{n_1^{0.25}}\right)\log_{10} \left(\frac{Re_v}{10000}\right) .. math:: F_{R,2} = \min(\frac{0.026}{F_L}\sqrt{n_1 Re_v},\; 1) .. math:: n_1 = \frac{N_2}{\left(\frac{C}{d^2}\right)^2} .. math:: F_R = F_{R,2} \text{ if Rev < 10 else } \min(F_{R,1a}, F_{R,2}) Otherwise : .. math:: F_{R,3a} = 1 + \left(\frac{0.33F_L^{0.5}}{n_2^{0.25}}\right)\log_{10} \left(\frac{Re_v}{10000}\right) .. math:: F_{R,4} = \frac{0.026}{F_L}\sqrt{n_2 Re_v} .. math:: n_2 = 1 + N_{32}\left(\frac{C}{d}\right)^{2/3} .. math:: F_R = F_{R,4} \text{ if Rev < 10 else } \min(F_{R,3a}, F_{R,4}) Parameters ---------- FL : float Liquid pressure recovery factor of a control valve without attached fittings [] C : float Metric Kv valve flow coefficient (flow rate of water at a pressure drop of 1 bar) [m^3/hr] d : float Diameter of the valve [m] Rev : float Valve reynolds number [-] full_trim : bool Whether or not the valve has full trim Returns ------- FR : float Reynolds number factor for laminar or transitional flow [] Examples -------- In Example 4, compressible flow with small flow trim sized for gas flow (Cv in the problem was converted to Kv here to make FR match with N32, N2): >>> Reynolds_factor(FL=0.98, C=0.015483, d=15., Rev=1202., full_trim=False) 0.7148753122302025 References ---------- .. [1] IEC 60534-2-1 / ISA-75.01.01-2007 ''' if full_trim: n1 = N2/(min(C/(d*d), 0.04))**2 # C/d**2 must not exceed 0.04 FR_1a = 1.0 + (0.33*sqrt(FL))/sqrt(sqrt(n1))*log10(Rev/10000.) FR_2 = 0.026/FL*sqrt(n1*Rev) if Rev < 10.0: FR = FR_2 else: FR = min(FR_2, FR_1a) else: n2 = 1 + N32*(C/d**2)**(2/3.) FR_3a = 1 + (0.33*sqrt(FL))/sqrt(sqrt(n2))*log10(Rev/10000.) FR_4 = min(0.026/FL*sqrt(n2*Rev), 1) if Rev < 10: FR = FR_4 else: FR = min(FR_3a, FR_4) return FR def size_control_valve_l(rho, Psat, Pc, mu, P1, P2, Q, D1=None, D2=None, d=None, FL=0.9, Fd=1, allow_choked=True, allow_laminar=True, full_output=False): r'''Calculates flow coefficient of a control valve passing a liquid according to IEC 60534. Uses a large number of inputs in SI units. Note the return value is not standard SI. All parameters are required. This sizing model does not officially apply to liquid mixtures, slurries, non-Newtonian fluids, or liquid-solid conveyance systems. For details of the calculations, consult [1]_. Parameters ---------- rho : float Density of the liquid at the inlet [kg/m^3] Psat : float Saturation pressure of the fluid at inlet temperature [Pa] Pc : float Critical pressure of the fluid [Pa] mu : float Viscosity of the fluid [Pa*s] P1 : float Inlet pressure of the fluid before valves and reducers [Pa] P2 : float Outlet pressure of the fluid after valves and reducers [Pa] Q : float Volumetric flow rate of the fluid [m^3/s] D1 : float, optional Diameter of the pipe before the valve [m] D2 : float, optional Diameter of the pipe after the valve [m] d : float, optional Diameter of the valve [m] FL : float, optional Liquid pressure recovery factor of a control valve without attached fittings (normally 0.8-0.9 at full open and decreasing as opened further to below 0.5; use default very cautiously!) [] Fd : float, optional Valve style modifier (0.1 to 1; varies tremendously depending on the type of valve and position; do not use the default at all!) [] allow_choked : bool, optional Overrides the automatic transition into the choked regime if this is False and returns as if choked flow does not exist allow_laminar : bool, optional Overrides the automatic transition into the laminar regime if this is False and returns as if laminar flow does not exist full_output : bool, optional If True, returns intermediate calculation values as well as Kv in the form of a dictionary containing 'Kv', 'Rev', 'choked', 'FL', 'FLP', 'FR', 'FP', and 'laminar'. Some may be None if they are not used in the calculation. Returns ------- Kv : float Metric Kv valve flow coefficient (flow rate of water at a pressure drop of 1 bar) [m^3/hr] Notes ----- It is possible to use this model without any diameters specified; in that case, turbulent flow is assumed. Choked flow can still be modeled. This is not recommended. All three diameters need to be None for this to work. `FL` and `Fd` are not used by the models when the diameters are not specified. Examples -------- From [1]_, matching example 1 for a globe, parabolic plug, flow-to-open valve. >>> size_control_valve_l(rho=965.4, Psat=70.1E3, Pc=22120E3, mu=3.1472E-4, ... P1=680E3, P2=220E3, Q=0.1, D1=0.15, D2=0.15, d=0.15, ... FL=0.9, Fd=0.46) 164.9954763704956 From [1]_, matching example 2 for a ball, segmented ball, flow-to-open valve. >>> size_control_valve_l(rho=965.4, Psat=70.1E3, Pc=22120E3, mu=3.1472E-4, ... P1=680E3, P2=220E3, Q=0.1, D1=0.1, D2=0.1, d=0.1, ... FL=0.6, Fd=0.98) 238.05817216710483 References ---------- .. [1] IEC 60534-2-1 / ISA-75.01.01-2007 ''' if full_output: ans = {'FLP': None, 'FP': None, 'FR': None} # Pa to kPa, according to constants in standard P1, P2, Psat, Pc = P1/1000., P2/1000., Psat/1000., Pc/1000. Q = Q*3600. # m^3/s to m^3/hr, according to constants in standard nu = mu/rho # kinematic viscosity used in standard MAX_C_POSSIBLE = 1E40 # Quit iterations if C reaches this high dP = P1 - P2 FF = FF_critical_pressure_ratio_l(Psat=Psat, Pc=Pc) choked = is_choked_turbulent_l(dP=dP, P1=P1, Psat=Psat, FF=FF, FL=FL) if choked and allow_choked: # Choked flow, equation 3 C = Q/N1/FL*sqrt(rho/rho0/(P1 - FF*Psat)) else: # non-choked flow, eq 1 C = Q/N1*sqrt(rho/rho0/dP) if D1 is None and D2 is None and d is None: # Assume turbulent if no diameters are provided, no other calculations Rev = 1e5 else: # m to mm, according to constants in standard D1, D2, d = D1*1000., D2*1000., d*1000. Rev = Reynolds_valve(nu=nu, Q=Q, D1=D1, FL=FL, Fd=Fd, C=C) # normal calculation path if (Rev > 10000 or not allow_laminar) and (D1 != d or D2 != d): # liquid, using Fp and FLP FP = 1.0 Ci = C MAX_ITER = 20 def iterate_piping_turbulent_l(Ci, iterations): loss = loss_coefficient_piping(d, D1, D2) FP = 1.0/sqrt(1 + loss/N2*(Ci/d**2)**2) if d > D1: loss_upstream = 0.0 else: loss_upstream = loss_coefficient_piping(d, D1) FLP = FL*1.0/sqrt(1 + FL**2/N2*loss_upstream*(Ci/d**2)**2) choked = is_choked_turbulent_l(dP, P1, Psat, FF, FLP=FLP, FP=FP) if choked: # Choked flow with piping, equation 4 C = Q/N1/FLP*sqrt(rho/rho0/(P1-FF*Psat)) else: # Non-Choked flow with piping, equation 5 C = Q/N1/FP*sqrt(rho/rho0/dP) if Ci/C < 0.99 and iterations < MAX_ITER and Ci < MAX_C_POSSIBLE: C = iterate_piping_turbulent_l(C, iterations+1) if MAX_ITER == iterations or Ci >= MAX_C_POSSIBLE: ans['warning'] = 'Not converged in inner loop' if full_output: ans['FLP'] = FLP ans['FP'] = FP return C C = iterate_piping_turbulent_l(Ci, 0) elif Rev <= 10000 and allow_laminar: # Laminar def iterate_piping_laminar_l(C): Ci = 1.3*C Rev = Reynolds_valve(nu=nu, Q=Q, D1=D1, FL=FL, Fd=Fd, C=Ci) if Ci/d**2 > 0.016*N18: FR = Reynolds_factor(FL=FL, C=Ci, d=d, Rev=Rev, full_trim=False) else: FR = Reynolds_factor(FL=FL, C=Ci, d=d, Rev=Rev, full_trim=True) if C/FR >= Ci: Ci = iterate_piping_laminar_l(Ci) # pragma: no cover if full_output: ans['Rev'] = Rev ans['FR'] = FR return Ci C = iterate_piping_laminar_l(C) if full_output: ans['FF'] = FF ans['choked'] = choked ans['Kv'] = C ans['laminar'] = Rev <= 10000 # For the laminar case this is already set and needs to not be overwritten if 'Rev' not in ans: ans['Rev'] = Rev return ans else: # return C, choked, laminar, FF, FR, Rev, FP, FLP, warning return C def size_control_valve_g(T, MW, mu, gamma, Z, P1, P2, Q, D1=None, D2=None, d=None, FL=0.9, Fd=1, xT=0.7, allow_choked=True, allow_laminar=True, full_output=False): r'''Calculates flow coefficient of a control valve passing a gas according to IEC 60534. Uses a large number of inputs in SI units. Note the return value is not standard SI. All parameters are required. For details of the calculations, consult [1]_. Note the inlet gas flow conditions. Parameters ---------- T : float Temperature of the gas at the inlet [K] MW : float Molecular weight of the gas [g/mol] mu : float Viscosity of the fluid at inlet conditions [Pa*s] gamma : float Specific heat capacity ratio [-] Z : float Compressibility factor at inlet conditions, [-] P1 : float Inlet pressure of the gas before valves and reducers [Pa] P2 : float Outlet pressure of the gas after valves and reducers [Pa] Q : float Volumetric flow rate of the gas at *273.15 K* and 1 atm specifically [m^3/s] D1 : float, optional Diameter of the pipe before the valve [m] D2 : float, optional Diameter of the pipe after the valve [m] d : float, optional Diameter of the valve [m] FL : float, optional Liquid pressure recovery factor of a control valve without attached fittings (normally 0.8-0.9 at full open and decreasing as opened further to below 0.5; use default very cautiously!) [] Fd : float, optional Valve style modifier (0.1 to 1; varies tremendously depending on the type of valve and position; do not use the default at all!) [] xT : float, optional Pressure difference ratio factor of a valve without fittings at choked flow (increasing to 0.9 or higher as the valve is closed further and decreasing to 0.1 or lower as the valve is opened further; use default very cautiously!) [-] allow_choked : bool, optional Overrides the automatic transition into the choked regime if this is False and returns as if choked flow does not exist allow_laminar : bool, optional Overrides the automatic transition into the laminar regime if this is False and returns as if laminar flow does not exist full_output : bool, optional If True, returns intermediate calculation values as well as Kv in the form of a dictionary containing 'Kv', 'Rev', 'choked', 'Y', 'FR', 'FP', 'xTP', and 'laminar'. Some may be None if they are not used in the calculation. Returns ------- Kv : float Metric Kv valve flow coefficient (flow rate of water at a pressure drop of 1 bar) [m^3/hr] Notes ----- It is possible to use this model without any diameters specified; in that case, turbulent flow is assumed. Choked flow can still be modeled. This is not recommended. All three diameters need to be None for this to work. `FL` and `Fd` are not used by the models when the diameters are not specified, but `xT` definitely is used by the model. When this model does not converge, the result is normally because of the specified delta P being less than that caused by the piping diameter changes. Examples -------- From [1]_, matching example 3 for non-choked gas flow with attached fittings and a rotary, eccentric plug, flow-to-open control valve: >>> size_control_valve_g(T=433., MW=44.01, mu=1.4665E-4, gamma=1.30, ... Z=0.988, P1=680E3, P2=310E3, Q=38/36., D1=0.08, D2=0.1, d=0.05, ... FL=0.85, Fd=0.42, xT=0.60) 72.5866454539105 From [1]_, roughly matching example 4 for a small flow trim sized tapered needle plug valve. Difference is 3% and explained by the difference in algorithms used. >>> size_control_valve_g(T=320., MW=39.95, mu=5.625E-5, gamma=1.67, Z=1.0, ... P1=2.8E5, P2=1.3E5, Q=0.46/3600., D1=0.015, D2=0.015, d=0.015, FL=0.98, ... Fd=0.07, xT=0.8) 0.016498765335995726 References ---------- .. [1] IEC 60534-2-1 / ISA-75.01.01-2007 ''' MAX_C_POSSIBLE = 1E40 # Quit iterations if C reaches this high # Pa to kPa, according to constants in standard P1, P2 = P1*1e-3, P2*1e-3 Q = Q*3600. # m^3/s to m^3/hr, according to constants in standard # Convert dynamic viscosity to kinematic viscosity Vm = Z*R*T/(P1*1000) rho = MW*1e-3/Vm nu = mu/rho # kinematic viscosity used in standard dP = P1 - P2 Fgamma = gamma/1.40 x = dP/P1 Y = max(1 - x/(3*Fgamma*xT), 2/3.) choked = is_choked_turbulent_g(x, Fgamma, xT) if choked and allow_choked: # Choked, and flow coefficient from eq 14a C = Q/(N9*P1*Y)*sqrt(MW*T*Z/xT/Fgamma) else: # Non-choked, and flow coefficient from eq 8a C = Q/(N9*P1*Y)*sqrt(MW*T*Z/x) if full_output: # numba: delete ans = {'FP': None, 'xTP': None, 'FR': None, 'choked': choked, 'Y': Y} # numba: delete if D1 is None and D2 is None and d is None: # Assume turbulent if no diameters are provided, no other calculations Rev = 1e5 if full_output: # numba: delete ans['Rev'] = None # numba: delete else: # m to mm, according to constants in standard D1, D2, d = D1*1000., D2*1000., d*1000. # Convert diameters to mm which is used in the standard Rev = Reynolds_valve(nu=nu, Q=Q, D1=D1, FL=FL, Fd=Fd, C=C) if full_output: # numba: delete ans['Rev'] = Rev # numba: delete if (Rev > 10000 or not allow_laminar) and (D1 != d or D2 != d): # gas, using xTP and FLP FP = 1. MAX_ITER = 20 def iterate_piping_coef_g(Ci, iterations): loss = loss_coefficient_piping(d, D1, D2) FP = 1.0/sqrt(1. + loss/N2*(Ci/d**2)**2) loss_upstream = loss_coefficient_piping(d, D1) xTP = xT/FP**2/(1 + xT*loss_upstream/N5*(Ci/d**2)**2) choked = is_choked_turbulent_g(x, Fgamma, xTP=xTP) if choked: # Choked flow with piping, equation 17a C = Q/(N9*FP*P1*Y)*sqrt(MW*T*Z/xTP/Fgamma) else: # Non-choked flow with piping, equation 11a C = Q/(N9*FP*P1*Y)*sqrt(MW*T*Z/x) if Ci/C < 0.99 and iterations < MAX_ITER and Ci < MAX_C_POSSIBLE: C = iterate_piping_coef_g(C, iterations+1) if full_output: # numba: delete ans['xTP'] = xTP # numba: delete ans['FP'] = FP # numba: delete ans['choked'] = choked # numba: delete if MAX_ITER == iterations or Ci >= MAX_C_POSSIBLE: # numba: delete ans['warning'] = 'Not converged in inner loop' # numba: delete return C # def err_piping_coeff(Ci): # loss = loss_coefficient_piping(d, D1, D2) # FP = (1. + loss/N2*(Ci/d**2)**2)**-0.5 # loss_upstream = loss_coefficient_piping(d, D1) # xTP = xT/FP**2/(1 + xT*loss_upstream/N5*(Ci/d**2)**2) # choked = is_choked_turbulent_g(x, Fgamma, xTP=xTP) # if choked: # # Choked flow with piping, equation 17a # C = Q/(N9*FP*P1*Y)*(MW*T*Z/xTP/Fgamma)**0.5 # else: # # Non-choked flow with piping, equation 11a # C = Q/(N9*FP*P1*Y)*(MW*T*Z/x)**0.5 # return C - Ci # import matplotlib.pyplot as plt # from fluids.numerics import linspace # Cs = linspace(C/50, C*50, 5000) # errs = [err_piping_coeff(C_test) for C_test in Cs] # plt.plot(Cs, errs) # plt.show() C = iterate_piping_coef_g(C, 0) elif Rev <= 10000 and allow_laminar: # Laminar; def iterate_piping_laminar_g(C): Ci = 1.3*C Rev = Reynolds_valve(nu=nu, Q=Q, D1=D1, FL=FL, Fd=Fd, C=Ci) if Ci/d**2 > 0.016*N18: FR = Reynolds_factor(FL=FL, C=Ci, d=d, Rev=Rev, full_trim=False) else: FR = Reynolds_factor(FL=FL, C=Ci, d=d, Rev=Rev, full_trim=True) if C/FR >= Ci: Ci = iterate_piping_laminar_g(Ci) if full_output: # numba: delete ans['FR'] = FR # numba: delete ans['Rev'] = Rev # numba: delete return Ci C = iterate_piping_laminar_g(C) if full_output: # numba: delete ans['Kv'] = C # numba: delete ans['laminar'] = Rev <= 10000 # numba: delete ans['choked'] = choked # numba: delete return ans # numba: delete return C # Valve data from Emerson Valve Handbook 5E # Quick opening valve data, spline fit, and interpolating function opening_quick = [0.0, 0.0136, 0.02184, 0.03256, 0.04575, 0.06221, 0.07459, 0.0878, 0.10757, 0.12654, 0.14301, 0.16032, 0.18009, 0.18999, 0.20233, 0.23105, 0.25483, 0.28925, 0.32365, 0.36541, 0.42188, 0.46608, 0.53319, 0.61501, 0.7034, 0.78033, 0.84415, 0.91944, 1.000] frac_CV_quick = [0.0, 0.04984, 0.07582, 0.12044, 0.16614, 0.21707, 0.26998, 0.32808, 0.39353, 0.46516, 0.52125, 0.58356, 0.64798, 0.68845, 0.72277, 0.76565, 0.79399, 0.82459, 0.84589, 0.86732, 0.88078, 0.89399, 0.90867, 0.92053, 0.93973, 0.95872, 0.96817, 0.98611, 1.0] opening_quick_tck = implementation_optimize_tck([[0.0, 0.0, 0.0, 0.0, 0.02184, 0.03256, 0.04575, 0.06221, 0.07459, 0.0878, 0.10757, 0.12654, 0.14301, 0.16032, 0.18009, 0.18999, 0.20233, 0.23105, 0.25483, 0.28925, 0.32365, 0.36541, 0.42188, 0.46608, 0.53319, 0.61501, 0.7034, 0.78033, 0.84415, 1.0, 1.0, 1.0, 1.0], [-3.2479258181113327e-19, 0.037650956835178835, 0.054616164261637117, 0.12657862552611354, 0.17115105822542115, 0.2075233903194021, 0.27084055195333684, 0.34208963001568016, 0.38730839943796663, 0.4656002247400036, 0.5196995880922897, 0.5907033063634928, 0.6304293931726886, 0.6953064258075168, 0.7382935002453699, 0.7631579537132379, 0.7997961180795559, 0.8262370617883222, 0.8471954722933543, 0.873096858463145, 0.8776128736976467, 0.897647305294458, 0.9105672165523071, 0.9192771703370824, 0.9377349743236904, 0.9603716623033031, 0.9688863605959851, 0.9980062718267431, 1.0, 0.0, 0.0, 0.0, 0.0], 3]) Cv_char_quick_opening = lambda opening: float(splev(opening, opening_quick_tck)) opening_linear = [0., 1.0] frac_CV_linear = [0, 1] Cv_char_linear = lambda opening: interp(opening, opening_linear, frac_CV_linear) # Equal opening valve data, spline fit, and interpolating function opening_equal = [0.0, 0.05523, 0.09287, 0.15341, 0.18942, 0.22379, 0.25816, 0.29582, 0.33348, 0.34985, 0.3826, 0.45794, 0.49235, 0.51365, 0.54479, 0.57594, 0.60218, 0.62843, 0.77628, 0.796, 0.83298, 0.86995, 0.90936, 0.95368, 1.00] frac_CV_equal = [0.0, 0.00845, 0.01339, 0.01877, 0.02579, 0.0349, 0.04189, 0.05528, 0.07079, 0.07533, 0.09074, 0.13444, 0.15833, 0.17353, 0.20159, 0.23388, 0.26819, 0.30461, 0.60113, 0.64588, 0.72583, 0.80788, 0.87519, 0.94999, 1.] opening_equal_tck = implementation_optimize_tck([[0.0, 0.0, 0.0, 0.0, 0.09287, 0.15341, 0.18942, 0.22379, 0.25816, 0.29582, 0.33348, 0.34985, 0.3826, 0.45794, 0.49235, 0.51365, 0.54479, 0.57594, 0.60218, 0.62843, 0.77628, 0.796, 0.83298, 0.86995, 0.90936, 1.0, 1.0, 1.0, 1.0], [1.3522591106779132e-19, 0.004087873896711868, 0.014374150571122216, 0.016455484312674015, 0.024946845435605228, 0.03592972456181881, 0.040710119644626126, 0.054518468768197687, 0.06976905178508139, 0.07587146190282387, 0.0985485829020452, 0.1238160142641967, 0.15558350087382017, 0.17487348629353283, 0.20157507610951217, 0.22995771158118564, 0.2683886931491415, 0.3574766835730407, 0.5027678906008036, 0.659729970241158, 0.7233389559355903, 0.8155475382785987, 0.8983628328699896, 0.9871204658597236, 1.0, 0.0, 0.0, 0.0, 0.0], 3]) Cv_char_equal_percentage = lambda opening: float(splev(opening, opening_equal_tck)) def convert_flow_coefficient(flow_coefficient, old_scale, new_scale): """Convert from one flow coefficient scale to another; supports the `Kv` `Cv`, and `Av` scales. Other scales are `Qn` and `Cg`, but clear definitions have yet to be found. Parameters ---------- flow_coefficient : float Value of the flow coefficient to be converted, expressed in the original scale. old_scale : str String specifying the original scale; one of 'Av', 'Cv', or 'Kv', [-] new_scale : str String specifying the new scale; one of 'Av', 'Cv', or 'Kv', [-] Returns ------- converted_flow_coefficient : float Flow coefficient converted to the specified scale. Notes ----- `Qn` is a scale based on a flow of air in units of L/minute as air travels through a valve and loses one bar of pressure (initially 7 bar absolute, to 6 bar absolute). No consistent conversion factors have been found and those from theory do not match what have been found. Some uses of `Qn` use its flow rate as in normal (STP reference conditions) flow rate of air; others use something like the 7 bar absolute condition. Examples -------- >>> convert_flow_coefficient(10, 'Kv', 'Av') 0.0002776532068951358 """ # Convert from `old_scale` to Kv if old_scale == 'Cv': Kv = Cv_to_Kv(flow_coefficient) elif old_scale == 'Kv': Kv = flow_coefficient elif old_scale == 'Av': Cv = flow_coefficient/(sqrt(rho0/psi)*gallon/minute) Kv = Cv_to_Kv(Cv) else: raise NotImplementedError("Supported scales are 'Cv', 'Kv', and 'Av'") if new_scale == 'Cv': ans = Kv_to_Cv(Kv) elif new_scale == 'Kv': ans = Kv elif new_scale == 'Av': Cv = Kv_to_Cv(Kv) ans = Cv*(sqrt(rho0/psi)*gallon/minute) else: raise NotImplementedError("Supported scales are 'Cv', 'Kv', and 'Av'") return ans # Third octave center frequency fi Hz fis_l_2015 = [12.5, 16.0, 20.0, 25.0, 31.5, 40.0, 50.0, 63.0, 80.0, 100.0, 125.0, 160.0, 200.0, 250.0, 315.0, 400.0, 500.0, 630.0, 800.0, 1000.0, 1250.0, 1600.0, 2000.0, 2500.0, 3150.0, 4000.0, 5000.0, 6300.0, 8000.0, 10000.0, 12500.0, 16000.0, 20000.0] #fis_l_2015_inv = [1.0/fi for fi in fis_l_2015] #fis_l_2015_1_5 = [fi**1.5 for fi in fis_l_2015] #fis_l_2015_n1_5 = [fi**-1.5 for fi in fis_l_2015] fis_l_2015_inv = [0.08, 0.0625, 0.049999999999999996, 0.04000000000000001, 0.031746031746031744, 0.025, 0.02, 0.01587301587301587, 0.012499999999999999, 0.010000000000000002, 0.008, 0.00625, 0.005, 0.003999999999999999, 0.003174603174603174, 0.0025000000000000005, 0.002, 0.0015873015873015873, 0.00125, 0.0009999999999999998, 0.0008, 0.0006250000000000001, 0.0005, 0.0004, 0.0003174603174603174, 0.00024999999999999995, 0.0002, 0.00015873015873015873, 0.000125, 0.0001, 8e-05, 6.249999999999999e-05, 5e-05] fis_l_2015_1_5 = [44.19417382415922, 64.0, 89.44271909999159, 125.0, 176.79331152506873, 252.98221281347034, 353.5533905932738, 500.04699779120756, 715.5417527999327, 1000.0, 1397.5424859373686, 2023.8577025077627, 2828.42712474619, 3952.847075210474, 5590.695395029137, 8000.0, 11180.339887498949, 15812.874501494027, 22627.41699796952, 31622.776601683792, 44194.17382415922, 64000.0, 89442.71909999159, 125000.0, 176793.3115250687, 252982.21281347034, 353553.39059327374, 500046.9977912076, 715541.7527999327, 1000000.0, 1397542.4859373686, 2023857.7025077627, 2828427.12474619] fis_l_2015_n1_5 = [0.02262741699796952, 0.015625, 0.011180339887498947, 0.008000000000000002, 0.00565632258015713, 0.003952847075210475, 0.00282842712474619, 0.001999812026503847, 0.0013975424859373684, 0.0010000000000000002, 0.0007155417527999327, 0.0004941058844013093, 0.00035355339059327376, 0.0002529822128134703, 0.00017886862533936855, 0.00012500000000000003, 8.944271909999159e-05, 6.323960895949173e-05, 4.419417382415922e-05, 3.162277660168379e-05, 2.2627416997969522e-05, 1.5625000000000004e-05, 1.1180339887498949e-05, 8.000000000000001e-06, 5.6563225801571285e-06, 3.9528470752104736e-06, 2.8284271247461903e-06, 1.9998120265038475e-06, 1.3975424859373686e-06, 1.0000000000000002e-06, 7.155417527999328e-07, 4.941058844013092e-07, 3.535533905932738e-07] fis_l_2015_3 = [1953.125, 4096.0, 8000.0, 15625.0, 31255.875, 64000.0, 125000.0, 250047.0, 512000.0, 1000000.0, 1953125.0, 4096000.0, 8000000.0, 15625000.0, 31255875.0, 64000000.0, 125000000.0, 250047000.0, 512000000.0, 1000000000.0, 1953125000.0, 4096000000.0, 8000000000.0, 15625000000.0, 31255875000.0, 64000000000.0, 125000000000.0, 250047000000.0, 512000000000.0, 1000000000000.0, 1953125000000.0, 4096000000000.0, 8000000000000.0] fis_l_2015_17 = [73.2397784872531, 111.43047210190387, 162.83621261476173, 237.95674233948478, 352.4746934040807, 529.0564156396547, 773.1237367774792, 1145.1936574895758, 1718.9093656438004, 2511.88643150958, 3670.6841971500585, 5584.753005453414, 8161.143093473476, 11926.088141608398, 17665.581651081215, 26515.632138719888, 38747.97468870842, 57395.64411646984, 86149.54298230256, 125892.54117941669, 183970.00582889825, 279900.6909294791, 409026.07302542904, 597720.3123687729, 885376.3998122095, 1328929.6319483411, 1941999.0242893337, 2876596.4096988947, 4317705.1125554005, 6309573.44480193, 9220341.829177868, 14028265.297730776, 20499864.602104142] #fis_l_2015_inv, fis_l_2015_1_5, fis_l_2015_17, fis_l_2015_n1_5, fis_l_2015_3 = [], [], [], [], [] #for fi in fis_l_2015: # fi_rt_inv = 1.0/sqrt(fi) # fis_l_2015_inv.append(fi_rt_inv*fi_rt_inv) # fis_l_2015_1_5.append(fi*fi*fi_rt_inv) # fis_l_2015_n1_5.append(fi_rt_inv*fi_rt_inv*fi_rt_inv) # fis_l_2015_3.append(fi*fi*fi) # fis_l_2015_17.append(fi**1.7) fis_length = 33 # dLa(fi), dB A_weights_l_2015 = [-63.4, -56.7, -50.5, -44.7, -39.4, -34.6, -30.2, -26.2, -22.5, -19.1, -16.1, -13.4, -10.9, -8.6, -6.6, -4.8, -3.2, -1.9, -0.8, 0.0, 0.6, 1.0, 1.2, 1.3, 1.2, 1.0, 0.5, -0.1, -1.1, -2.5, -4.3, -6.6, -9.3] def control_valve_noise_l_2015(m, P1, P2, Psat, rho, c, Kv, d, Di, FL, Fd, t_pipe, rho_pipe=7800.0, c_pipe=5000.0, rho_air=1.2, c_air=343.0, xFz=None, An=-4.6): r'''Calculates the sound made by a liquid flowing through a control valve according to the standard IEC 60534-8-4 (2015) [1]_. Parameters ---------- m : float Mass flow rate of liquid through the control valve, [kg/s] P1 : float Inlet pressure of the fluid before valves and reducers [Pa] P2 : float Outlet pressure of the fluid after valves and reducers [Pa] Psat : float Saturation pressure of the fluid at inlet temperature [Pa] rho : float Density of the liquid at the inlet [kg/m^3] c : float Speed of sound of the liquid at the inlet conditions [m/s] Kv : float Metric Kv valve flow coefficient (flow rate of water at a pressure drop of 1 bar) [m^3/hr] d : float Diameter of the valve [m] Di : float Internal diameter of the pipe before and after the valve [m] FL : float, optional Liquid pressure recovery factor of a control valve without attached fittings (normally 0.8-0.9 at full open and decreasing as opened further to below 0.5) [-] Fd : float, optional Valve style modifier [-] t_pipe : float Wall thickness of the pipe after the valve, [m] rho_pipe : float, optional Density of the pipe wall material at flowing conditions, [kg/m^3] c_pipe : float, optional Speed of sound of the pipe wall material at flowing conditions, [m/s] rho_air : float, optional Density of the air surrounding the valve and pipe wall, [kg/m^3] c_air : float, optional Speed of sound of the air surrounding the valve and pipe wall, [m/s] xFz : float, optional If specified, this value `xFz` is used instead of estimated; the calculation is sensitive to this value, [-] An : float, optional Valve correction factor for acoustic efficiency Returns ------- LpAe1m : float A weighted sound pressure level 1 m from the pipe wall, 1 m distance dowstream of the valve (at reference sound pressure level 2E-5), [dB] Notes ----- For formulas see [1]_. This takes on the order of 100 us to compute. This model can also tell if noise is being produced in a valve just due to turbulent flow, or cavitation. For values of `An`, see [1]_; it is normally -4.6 for global valves, -4.3 for butterfly valves, and -4.0 for expanders. This model was checked against three examples in [1]_; they match to all given decimals. A formula is given in [1]_ for multihole trim valves to estimate `xFz` as well; this is not implemented here and `xFz` must be calculated by the user separately. The formula is .. math:: x_{Fz} = \left(4.5 + 1650\frac{N_0d_H^2}{F_L}\right)^{-1/2} Where `N0` is the number of open channels and `dH` is the multihole trim hole diameter. Examples -------- >>> control_valve_noise_l_2015(m=40, P1=1E6, P2=6.5E5, Psat=2.32E3, ... rho=997, c=1400, Kv=77.848, d=0.1, Di=0.1071, FL=0.92, Fd=0.42, ... t_pipe=0.0036, rho_pipe=7800.0, c_pipe=5000.0,rho_air=1.293, ... c_air=343.0, An=-4.6) 81.58200097996 References ---------- .. [1] IEC 60534-8-4 : Industrial-Process Control Valves - Part 8-4: Noise Considerations - Prediction of Noise Generated by Hydrodynamic Flow. (2015) ''' # Convert Kv to Cv as C N34 = 1.17 # for Cv - conversion constant but not to many decimals N14 = 0.0046 C = Kv_to_Cv(Kv) xF = (P1-P2)/(P1-Psat) dPc = min(P1-P2, FL*FL*(P1 - Psat)) if xFz is None: xFz = 0.9*1.0/sqrt(1.0 + 3.0*Fd*sqrt(C/(N34*FL))) xFzp1 = xFz*sqrt(sqrt(sqrt(6E5/P1))) Dj = N14*Fd*sqrt(C*FL) Uvc = sqrt(2.0*dPc/rho)/FL Wm = 0.5*m*Uvc*Uvc*FL*FL cavitating = xF > xFzp1 eta_turb = 10.0**An*Uvc/c x0 = xF - xFzp1 x1 = xF/xFzp1 x2 = x1*x1 x1 = x2*x2*x1 if cavitating: eta_cav = 0.32*eta_turb*sqrt((P1 - P2)/(dPc*xFzp1))*exp(5.0*xFzp1)*sqrt((1.0 - xFzp1)/(1.0 - xF))*(x1)*x0*sqrt(x0) Wa = (eta_turb+eta_cav)*Wm else: Wa = eta_turb*Wm Lpi = 10.0*log10(3.2E9*Wa*rho*c/(Di*Di)) Stp = 0.036*FL*FL*C*Fd**0.75/(N34*xFzp1*sqrt(xFzp1)*d*d)*(1.0/(P1 - Psat))**0.57 f_p_turb = Stp*Uvc/Dj if cavitating: x3 = ((1.0 - xF)/(1.0 - xFzp1)) x4 = xFzp1/xF f_p_cav = 6.0*f_p_turb*x3*x3*x4*x4*sqrt(x4) f_p_cav_inv = 1.0/f_p_cav f_p_cav_inv_1_5 = f_p_cav_inv*sqrt(f_p_cav_inv) f_p_cav_inv_1_5_1_4 = 0.25*f_p_cav_inv_1_5 f_p_cav_1_5 = 1.0/f_p_cav_inv_1_5 eta_denom = 1.0/(eta_turb + eta_cav) t1 = eta_turb*eta_denom t2 = eta_cav*eta_denom fr = c_pipe/(pi*Di) fr_inv = 1.0/fr TL_fr = -10.0 - 10.0*log10(c_pipe*rho_pipe*t_pipe/(c_air*rho_air*Di)) t3 = - 10.0*log10((Di + 2.0*t_pipe + 2.0)/(Di + 2.0*t_pipe)) # F_cavs = [] # F_turbs = [] # LPis = [] # TL_fis = [] # L_pe1m_fis = [] LpAe1m_sum = 0.0 f_p_turb_inv = 1.0/f_p_turb f_p_turb_inv3 = f_p_turb_inv*f_p_turb_inv*f_p_turb_inv fr_inv_1_5 = fr_inv*sqrt(fr_inv) a_factor = ln_10_inv for i in range(fis_length): # for fi, fi_inv, fi_1_5, fi_1_5_inv, A in zip(fis_l_2015, fis_l_2015_inv, fis_l_2015_1_5, fis_l_2015_n1_5, A_weights_l_2015): # fi_inv = 1.0/fi # fi_turb_ratio = fis_l_2015[i]*f_p_turb_inv # fi_turb_ratio = fi*f_p_turb_inv F_turb = -.8 - log(0.25*f_p_turb_inv3*fis_l_2015_3[i] + fis_l_2015_inv[i]*f_p_turb)*a_factor # F_turbs.append(F_turb) if cavitating: # fi_cav_ratio = fi_1_5*f_p_cav_inv_1_5# (fi*f_p_cav_inv)**1.5 # F_cav = -.9 - log10(f_p_cav_inv_1_5_1_4*fis_l_2015_1_5[i] + fis_l_2015_n1_5[i]*f_p_cav_1_5) # 1.0/fi_cav_ratio, fi_1_5_inv*f_p_cav_1_5 F_cav_fact = 0.12589254117941673/(f_p_cav_inv_1_5_1_4*fis_l_2015_1_5[i] + fis_l_2015_n1_5[i]*f_p_cav_1_5) # 0.1258925411794167310 = 10**(-0.9) # 4.3429448190325175*log(x) -> 10*log10(x) LPif = (Lpi + 4.3429448190325175*log(t1*exp(ln_10*F_turb) + t2*F_cav_fact)) # Should be able to save 1 power in the above function somehow, combine the tow terms in exponent else: LPif = Lpi + F_turb*10.0 # LPis.append(LPif) # -8.685889638065035 = -20*log10(x) TL_fi = TL_fr - 8.685889638065035*log(fr*fis_l_2015_inv[i] + fis_l_2015_1_5[i]*fr_inv_1_5) # (fi*fr_inv)**1.5 # TL_fis.append(TL_fi) L_pe1m_fi = LPif + TL_fi + t3 # L_pe1m_fis.append(L_pe1m_fi) LpAe1m_sum += exp(0.23025850929940458*(L_pe1m_fi + A_weights_l_2015[i])) LpAe1m = 4.3429448190325175*log(LpAe1m_sum) return LpAe1m def control_valve_noise_g_2011(m, P1, P2, T1, rho, gamma, MW, Kv, d, Di, t_pipe, Fd, FL, FLP=None, FP=None, rho_pipe=7800.0, c_pipe=5000.0, P_air=101325.0, rho_air=1.2, c_air=343.0, An=-3.8, Stp=0.2, T2=None, beta=0.93): r'''Calculates the sound made by a gas flowing through a control valve according to the standard IEC 60534-8-3 (2011) [1]_. Parameters ---------- m : float Mass flow rate of gas through the control valve, [kg/s] P1 : float Inlet pressure of the gas before valves and reducers [Pa] P2 : float Outlet pressure of the gas after valves and reducers [Pa] T1 : float Inlet gas temperature, [K] rho : float Density of the gas at the inlet [kg/m^3] gamma : float Specific heat capacity ratio [-] MW : float Molecular weight of the gas [g/mol] Kv : float Metric Kv valve flow coefficient (flow rate of water at a pressure drop of 1 bar) [m^3/hr] d : float Diameter of the valve [m] Di : float Internal diameter of the pipe before and after the valve [m] t_pipe : float Wall thickness of the pipe after the valve, [m] Fd : float Valve style modifier (0.1 to 1; varies tremendously depending on the type of valve and position; do not use the default at all!) [-] FL : float Liquid pressure recovery factor of a control valve without attached fittings (normally 0.8-0.9 at full open and decreasing as opened further to below 0.5; use default very cautiously!) [-] FLP : float, optional Combined liquid pressure recovery factor with piping geometry factor, for a control valve with attached fittings [-] FP : float, optional Piping geometry factor [-] rho_pipe : float, optional Density of the pipe wall material at flowing conditions, [kg/m^3] c_pipe : float, optional Speed of sound of the pipe wall material at flowing conditions, [m/s] P_air : float, optional Pressure of the air surrounding the valve and pipe wall, [Pa] rho_air : float, optional Density of the air surrounding the valve and pipe wall, [kg/m^3] c_air : float, optional Speed of sound of the air surrounding the valve and pipe wall, [m/s] An : float, optional Valve correction factor for acoustic efficiency, [-] Stp : float, optional Strouhal number at the peak `fp`; between 0.1 and 0.3 typically, [-] T2 : float, optional Outlet gas temperature; assumed `T1` if not provided (a PH flash should be used to obtain this if possible), [K] beta : float, optional Valve outlet / expander inlet contraction coefficient, [-] Returns ------- LpAe1m : float A-weighted sound pressure level 1 m from the pipe wall, 1 m distance dowstream of the valve (at reference sound pressure level 2E-5), [dB] Notes ----- For formulas see [1]_. This takes on the order of 100 us to compute. For values of `An`, see [1]_. This model was checked against six examples in [1]_; they match to all given decimals. Several additional formulas are given for multihole trim valves, control valves with two or more fixed area stages, and multipath, multistage trim valves. Examples -------- >>> control_valve_noise_g_2011(m=2.22, P1=1E6, P2=7.2E5, T1=450, rho=5.3, ... gamma=1.22, MW=19.8, Kv=77.85, d=0.1, Di=0.2031, FL=None, FLP=0.792, ... FP=0.98, Fd=0.296, t_pipe=0.008, rho_pipe=8000.0, c_pipe=5000.0, ... rho_air=1.293, c_air=343.0, An=-3.8, Stp=0.2) 91.67702674629604 References ---------- .. [1] IEC 60534-8-3 : Industrial-Process Control Valves - Part 8-3: Noise Considerations - Control Valve Aerodynamic Noise Prediction Method." ''' k = gamma # alias C = Kv_to_Cv(Kv) N14 = 4.6E-3 N16 = 4.89E4 fs = 1.0 # structural loss factor reference frequency, Hz P_air_std = 101325.0 if T2 is None: T2 = T1 x = (P1 - P2)/P1 # FLP/FP when fittings attached FL_term = FLP/FP if FP is not None else FL P_vc = P1*(1.0 - x/FL_term**2) x_vcc = 1.0 - (2.0/(k + 1.0))**(k/(k - 1.0)) # mostly matches xc = FL_term**2*x_vcc alpha = (1.0 - x_vcc)/(1.0 - xc) xB = 1.0 - 1.0/alpha*(1.0/k)**(k/(k - 1.0)) xCE = 1.0 - 1.0/(22.0*alpha) # Regime determination check - should be ordered or won't work # assert xc < x_vcc # assert x_vcc < xB # assert xB < xCE if x <= xc: regime = 1 elif xc < x <= x_vcc: regime = 2 elif x_vcc < x <= xB: regime = 3 elif xB < x <= xCE: regime = 4 else: regime = 5 # print('regime', regime) Dj = N14*Fd*sqrt(C*(FL_term)) Mj5 = sqrt(2.0/(k - 1.0)*( 22.0**((k-1.0)/k) - 1.0 )) if regime == 1: Mvc = sqrt((2.0/(k-1.0)) *((1.0 - x/FL_term**2)**((1.0 - k)/k) - 1.0)) # Not match elif regime in (2, 3, 4): Mj = sqrt((2.0/(k-1.0))*((1.0/(alpha*(1.0-x)))**((k - 1.0)/k) - 1.0)) # Not match Mj = min(Mj, Mj5) # elif regime == 5: # pass if regime == 1: Tvc = T1*(1.0 - x/(FL_term)**2)**((k - 1.0)/k) cvc = sqrt(k*P1/rho*(1 - x/(FL_term)**2)**((k-1.0)/k)) Wm = 0.5*m*(Mvc*cvc)**2 else: Tvcc = 2.0*T1/(k + 1.0) cvcc = sqrt(2.0*k*P1/(k+1.0)/rho) Wm = 0.5*m*cvcc*cvcc # print('Wm', Wm) if regime == 1: fp = Stp*Mvc*cvc/Dj elif regime in (2, 3): fp = Stp*Mj*cvcc/Dj elif regime == 4: fp = 1.4*Stp*cvcc/Dj/sqrt(Mj*Mj - 1.0) elif regime == 5: fp = 1.4*Stp*cvcc/Dj/sqrt(Mj5*Mj5 - 1.0) fp_inv = 1.0/fp # print('fp', fp) if regime == 1: eta = 10.0**An*FL_term**2*(Mvc)**3 elif regime == 2: eta = 10.0**An*x/x_vcc*Mj**(6.6*FL_term*FL_term) elif regime == 3: eta = 10.0**An*Mj**(6.6*FL_term*FL_term) elif regime == 4: eta = 0.5*10.0**An*Mj*Mj*(sqrt(2.0))**(6.6*FL_term*FL_term) elif regime == 5: eta = 0.5*10.0**An*Mj5*Mj5*(sqrt(2.0))**(6.6*FL_term*FL_term) # print('eta', eta) Wa = eta*Wm rho2 = rho*(P2/P1) # Speed of sound c2 = sqrt(k*R*T2/(MW/1000.)) Mo = 4.0*m/(pi*d*d*rho2*c2) M2 = 4.0*m/(pi*Di*Di*rho2*c2) # print('M2', M2) Lg = 16.0*log10(1.0/(1.0 - min(M2, 0.3))) # dB if M2 > 0.3: Up = 4.0*m/(pi*rho2*Di*Di) UR = Up*Di*Di/(beta*d*d) WmR = 0.5*m*UR*UR*( (1.0 - d*d/(Di*Di))**2 + 0.2) fpR = Stp*UR/d MR = UR/c2 # Value listed in appendix here is wrong, "based on another # earlier standard. Calculation thereon is wrong". Assumed # correct, matches spreadsheet to three decimals. eta_R = 10**An*MR**3 WaR = eta_R*WmR L_piR = 10.0*log10((3.2E9)*WaR*rho2*c2/(Di*Di)) + Lg # print('Up', Up) # print('UR', UR) # print('WmR', WmR) # print('fpR', fpR) # print('MR', MR) # print('eta_R', eta_R, eta_R/8.8E-4) # print('WaR', WaR) # print('L_piR', L_piR) L_pi = 10.0*log10((3.2E9)*Wa*rho2*c2/(Di*Di)) + Lg # print('L_pi', L_pi) fr = c_pipe/(pi*Di) fo = 0.25*fr*(c2/c_air) fg = sqrt(3)*c_air**2/(pi*t_pipe*c_pipe) if d > 0.15: dTL = 0.0 elif 0.05 <= d <= 0.15: dTL = -16660.0*d**3 + 6370.0*d**2 - 813.0*d + 35.8 else: dTL = 9.0 # print(dTL, 'dTL') P_air_ratio = P_air/P_air_std LpAe1m_sum = 0.0 # LPis = [] # LPIRs = [] # L_pe1m_fis = [] for fi, A_weight in zip(fis_l_2015, A_weights_l_2015): # This gets adjusted when Ma > 0.3 fi_turb_ratio = fi*fp_inv t1 = 1.0 + (0.5*fi_turb_ratio)**2.5 t2 = 1.0 + (0.5/fi_turb_ratio)**1.7 # Formula forgot to use log10, but log10 is needed for the numbers Lpif = L_pi - 8.0 - 10.0*log10(t1*t2) # print(Lpif, 'Lpif') # LPis.append(Lpif) if M2 > 0.3: fiR_turb_ratio = fi/fpR t1 = 1.0 + (0.5*fiR_turb_ratio)**2.5 t2 = 1.0 + (0.5/fiR_turb_ratio)**1.7 # Again, log10 is missing LpiRf = L_piR - 8.0 - 10.0*log10(t1*t2) # LPIRs.append(LpiRf) LpiSf = 10.0*log10( 10**(0.1*Lpif) + 10.0**(0.1*LpiRf) ) if fi < fo: Gx = (fo/fr)**(2.0/3.0)*(fi/fo)**4.0 if fo < fg: Gy = (fo/fg) else: Gy = 1.0 else: if fi < fr: Gx = sqrt(fi/fr) else: Gx = 1.0 if fi < fg: Gy = fi/fg else: Gy = 1.0 eta_s = sqrt(0.01/fi) # print('eta_s', eta_s) # up to eta_s is good den = (rho2*c2 + 2.0*pi*t_pipe*fi*rho_pipe*eta_s)/(415.0*Gy) + 1.0 TL_fi = 10.0*log10(8.25E-7*(c2/(t_pipe*fi))**2*Gx/den*P_air_ratio) - dTL # Formula forgot to use log10, but log10 is needed for the numbers if M2 > 0.3: term = LpiSf else: term = Lpif L_pe1m_fi = term + TL_fi - 10.0*log10((Di + 2.0*t_pipe + 2.0)/(Di + 2.0*t_pipe)) # L_pe1m_fis.append(L_pe1m_fi) # print(L_pe1m_fi) LpAe1m_sum += 10.0**(0.1*(L_pe1m_fi + A_weight)) LpAe1m = 10.0*log10(LpAe1m_sum) return LpAe1m