diff --git a/PDSim/flow/ValveModels.py b/PDSim/flow/ValveModels.py new file mode 100644 index 0000000..308243b --- /dev/null +++ b/PDSim/flow/ValveModels.py @@ -0,0 +1,513 @@ +from __future__ import division, print_function + +import numpy as np +from math import log,tan,sin,cos,pi + +from PDSim.misc.datatypes import arraym, empty_arraym +from scipy.integrate import ode + +# def public: +# OUTPUT_VELOCITY +# OUTPUT_MA + +class ValveModel(object): + """ + + .. plot:: + + import matplotlib.pyplot as plt + import numpy as np + + fig=plt.figure() + ax1=fig.add_subplot(111) + ax1.fill(np.r_[0.2,0.4,0.4,0.2,0.2],np.r_[0,0,1,1,0],'grey') + ax1.text(0.41,0.66,r'$\\leftarrow p_{low}A_{valve}$',size=20,ha='left',va='center') + ax1.text(0.41,0.33,r'$\\leftarrow k_{valve}x$',size=20,ha='left',va='center') + ax1.text(0.19,0.66,r'$p_{high}A_{valve}\\rightarrow$',size=20,ha='right',va='center') + ax1.text(0.19,0.33,r'$\\frac{1}{2}C_D\\rho (V-\dot x)^2 A_{valve}\\rightarrow$',size=20,ha='right',va='center') + ax1.set_xlim(0.1-1,0.1+1) + ax1.axis('equal') + ax1.axis('off') + ax1.set_title('Pressure-dominant Free-body-diagram') + plt.show() + + Pressure-dominant Regime + + .. math:: + + M_{valve}\ddot x_{valve}+k_{valve}x_{valve} = (p_{high}-p_{low}) A_{valve}+\\frac{1}{2}C_D\\rho (V-\dot x_{valve})^2A_{valve} + + Two variables are :math:`x_2=\dot x_{valve}` and :math:`x_1=x_{valve}` where :math:`\ddot x_{valve}=\\frac{d}{dt}[\dot x_{valve}]` or :math:`\ddot x_{valve}=\dot x_2`. Thus the system of derivatives is + + .. math:: + + \mathbf f_{valves}=\\frac{d}{dt}\left[ \\begin{array}{c} \\dot x_{valve} \\\\ x_{valve} \end{array} \\right]=\left[ \\begin{array}{c} \\frac{d}{dt}[\dot x_{valve}] \\\\ \\frac{d}{dt}[x_{valve}] \end{array} \\right] + + Thus the system of equations is given by + + .. math:: + + \dot x_1 = x_2 + + M_{valve}\dot x_2+k_{valve}x_1 = (p_{high}-p_{low}) A_{valve}+\\frac{1}{2}C_D\\rho (V-x_2)^2A_{valve} + + which yields the solution for the derivatives of :math:`x_1` and :math:`x_2` of + + .. math:: + + \dot x_1 = x_2 + + \dot x_2 = \dfrac{(p_{high}-p_{low}) A_{valve}+\\frac{1}{2}C_D\\rho (V-x_2)^2A_{valve}-k_{valve}x_1}{M_{valve}} + + .. plot:: + + import matplotlib.pyplot as plt + import numpy as np + + fig=plt.figure() + ax1=fig.add_subplot(111) + ax1.fill(np.r_[0.2,0.4,0.4,0.2,0.2],np.r_[0,0,1,1,0],'grey') + ax1.text(0.41,0.66,r'$\leftarrow p_{low}A_{valve}$',size=20,ha='left',va='center') + ax1.text(0.41,0.33,r'$\leftarrow k_{valve}x$',size=20,ha='left',va='center') + ax1.text(0.19,0.8,r'$p_{low}A_{valve} \\rightarrow$',size=20,ha='right',va='center') + ax1.text(0.19,0.5,r'$\\frac{1}{2}C_D\\rho (V-\dot x)^2 A_{valve} \\rightarrow$',size=20,ha='right',va='center') + ax1.text(0.19,0.2,r'$\\rho (V-\dot x)^2 A_{port} \\rightarrow$',size=20,ha='right',va='center') + ax1.set_xlim(0.1-1,0.1+1) + ax1.axis('equal') + ax1.axis('off') + ax1.set_title('Mass-flux-dominant Free-body-diagram') + plt.show() + + And if mass-flux-dominated, force balance given by + + .. math:: + + M_{valve}\ddot x_{valve}+k_{valve}x_{valve} = \\frac{1}{2}C_D\\rho (\mathbf V-\dot x_{valve})^2 A_{valve}+\\rho (\mathbf V-\dot x_{valve})^2 A_{port} + + Which yields the solution for the system of derivatives of + + .. math:: + + \dot x_1 = x_2 + + \dot x_2= \dfrac{\\frac{1}{2}C_D\\rho (\mathbf V-x_2)^2 A_{valve}+\\rho (\mathbf V-x_2)^2 A_{port}-k_{valve}x_1}{M_{valve}} + """ + + def __init__(self, d_valve, d_port, C_D, rho_valve, x_stopper,m_eff,k_valve,x_tr,key_up, key_down): + ''' + # this function will be run when run the valve model first time + ''' + self.xv = empty_arraym(2) + self.rho_valve = rho_valve + self.d_valve = d_valve + self.d_port = d_port + self.A_valve = pi*d_valve**2/4.0 + self.A_port = pi*d_port**2/4.0 + self.m_eff = m_eff + self.C_D = C_D + self.k_valve = k_valve + self.x_stopper = x_stopper + self.key_up = key_up + self.key_down = key_down + self.x_tr = x_tr + + def get_States(self, Core): + """ + Core is the main model core, it contains information that + is needed for the flow models + + self.CVs=ControlVolumeCollection() + self.exists_keys = [self.keys[i] for i in self.exists_indices] + self.exists_indices = [i for i in self.indices if self.CVs[i].exists] + CVs : list of control volumes + + self.Tubes = TubeCollection() + + """ + exists_keys=Core.CVs.exists_keys # exists : bool ``True`` if control volume exists, ``False`` otherwise + Tubes_Nodes=Core.Tubes.Nodes + + for key, Statevar in [(self.key_up,'State_up'),(self.key_down,'State_down')]: + + # Update the pointers to the states for the ends of the flow path + if key in exists_keys: + setattr(self,Statevar,Core.CVs[key].State) + elif key in Tubes_Nodes: + setattr(self,Statevar,Tubes_Nodes[key]) + + def set_xv(self,xv): + + self.xv = xv.copy() + # if self.xv[0] < -1e-15 and self.xv.get_index(1) < 1e-15: + if self.xv[0] < 0 and self.xv[1] < 1e-15: + self.xv[0] = 0.0 + self.xv[1] = 0.0 + #If it predicts a valve opening greater than max opening, just use the max opening + elif self.xv[0] > self.x_stopper and self.xv[1] > 0: + self.xv[0] = self.x_stopper + self.xv[1] = 0.0 + + def _pressure_dominant(self, f, x,xdot,rho,V, deltap): + """ + set_index() in cython becomes standard array assignment + """ + + f[0] = xdot + + if abs(V-xdot) > 0: + f[1] = ((V-xdot)/abs(V-xdot)*0.5*self.C_D*rho*V**2*self.A_valve+deltap*self.A_valve-self.k_valve*x)/(self.m_eff) #xdot + else: + f[1] = (deltap*self.A_valve-self.k_valve*x)/(self.m_eff) #d(xdot)dt + + return + + def _flux_dominant(self, f,x, xdot,rho,V): + + f[0] = xdot #dxdt + + if abs(V-xdot) > 0: + f[1] = ((V-xdot)/abs(V-xdot)*0.5*self.C_D*rho*V**2*self.A_valve+(V-xdot)/abs(V-xdot)*rho*(V-xdot)**2*self.A_port-self.k_valve*x)/(self.m_eff) #xdot + else: + f[1] = (-self.k_valve*x)/(self.m_eff) #d(xdot)dt + + return + + def A(self): + + if self.xv is None: + print('self.xv is None') + x = self.xv[0] + if x >= self.x_tr: + return self.A_port # flux domain + else: + return pi*x*self.d_valve # pressure domain + + def flow_velocity(self, State_up, State_down): + """ + For a given set of states, and a known valve lift, first + check whether it is within the valve lift range, and then + calculate the flow velocity + """ + A = self.A() # From def A to get the calculation area + # print ' A in flow_velocity is:',A # A is 0 all the time + x = self.xv[0] + + if A > 0: + if x > self.x_tr: + return IsentropicNozzle(A, State_up, State_down, 2) # use the default mdot as our output + else: + # return x/self.x_tr*IsentropicNozzle(A, State_up, State_down, OUTPUT_VELOCITY) + return x/self.x_tr*IsentropicNozzle(A, State_up, State_down, 2) # can output velocity in cpython with public + else: + return 0.0 + + def derivs(self, Core): + """ + Return the position and velocity as an arraym for the valve + + Parameters + ---------- + Core : :class:`PDSimCore ` instance + + Returns + ------- + out_array : :class:`arraym ` instance + + """ + f = empty_arraym(2) + out_array = empty_arraym(2) + x = self.xv[0] + xdot = self.xv[1] + + self.get_States(Core) + + rho = self.State_up.get_rho() + p_high = self.State_up.get_p() + p_low = self.State_down.get_p() + deltap = (p_high - p_low)*1000 + + if deltap > 0: + V = self.flow_velocity(self.State_up, self.State_down) + else: + V = -self.flow_velocity(self.State_down, self.State_up) + + if x <= self.x_tr: + self._pressure_dominant(f,x,xdot,rho,V,deltap) + else: + self._flux_dominant(f,x,xdot,rho,V) + + omega = Core.omega + + out_array[0] = f[0]/omega + out_array[1] = f[1]/omega + + #TODO: check here (3)!!! + if abs(x) < 1e-15 and xdot < -1e-12: + # print 'stationary valve' + out_array[0] = 0.0 + out_array[1] = 0.0 + + return out_array #[dxdtheta, d(xdot)_dtheta] + + #TODO: check __cdict__ + def __dict__(self): + + # print 'Im in dict!' + items=['A_port','A_valve','d_valve','h_valve','d_port','m_eff','C_D','a_valve','l_valve', + 'rho_valve','k_valve','x_stopper','key_up','key_down','xv','State_up', 'State_down', + 'E'] + + return {item:getattr(self,item) for item in items} + + def __repr__(self): + """ + Return a representation of Valve Model for outputting to screen + """ + # print 'Im in repr!' + s='' + for item in self.__dict__(): + s += item+' : '+str(getattr(self,item))+'\n' + return s + + def __reduce__(self): + # print 'Im in reduce!' + return rebuildValveModel,(self.__dict__(),) #TODO: check __cdict__ + + +class ValveModel_kvalve(object): + """ + Mathison M.M., PhD Thesis, Purdue University, 2011 + + kvalve = a*exp(b*xv) + c + + """ + def __init__(self, valve_type,d_valve, d_port, C_D,cheta, x_stopper,m_valve,x_tr,key_up, key_down): + ''' + # this function will be run when run the valve model first time + ''' + self.xv = empty_arraym(2) + + self.valve_type = valve_type + self.d_valve = d_valve + self.d_port = d_port + self.A_valve = pi*d_valve**2/4.0 + self.A_port = pi*d_port**2/4.0 + self.m_valve = m_valve + self.C_D = C_D + self.cheta = cheta + self.x_stopper = x_stopper + self.key_up = key_up + self.key_down = key_down + self.x_tr = x_tr + + self.Bv = pi*(self.d_port * 1.14)**2/ 4 #Effective area of the discharging port + + + def get_States(self, Core): + """ + Core is the main model core, it contains information that + is needed for the flow models + + self.CVs=ControlVolumeCollection() + self.exists_keys = [self.keys[i] for i in self.exists_indices] + self.exists_indices = [i for i in self.indices if self.CVs[i].exists] + CVs : list of control volumes + + self.Tubes = TubeCollection() + + """ + exists_keys=Core.CVs.exists_keys # exists : bool ``True`` if control volume exists, ``False`` otherwise + Tubes_Nodes=Core.Tubes.Nodes + + + for key, Statevar in [(self.key_up,'State_up'),(self.key_down,'State_down')]: + + # Update the pointers to the states for the ends of the flow path + if key in exists_keys: + setattr(self,Statevar,Core.CVs[key].State) + elif key in Tubes_Nodes: + setattr(self,Statevar,Tubes_Nodes[key]) + + def set_xv(self,xv): + + self.xv = xv.copy() + + # if self.xv[0] < -1e-15 and self.xv.get_index(1) < 1e-15: + if self.xv[0] <= 0:# and self.xv[1] < 1e-15: + self.xv[0] = 0.0 + self.xv[1] = 0.0 + #If it predicts a valve opening greater than max opening, just use the max opening + elif self.xv[0] > self.x_stopper: # and self.xv[1] > 0: + self.xv[0] = self.x_stopper + self.xv[1] = 0.0 + + def _pressure_dominant(self, f, x,xdot,Dl,deltap): + """ + set_index() in cython becomes standard array assignment + """ + f[0] = xdot # Velovity of a valve + f[1] = Dl - 2*self.cheta * self.omega_n * xdot - self.omega_n**2*x #Acceleration of a valve + + return + + + def A(self): + + if self.xv is None: + print('self.xv is None') + x = self.xv[0] + if x <= 0.0: + return 0.0 + elif x >= self.x_tr: + return self.A_port # flux domain + else: + Av = self.Bv / np.sqrt(1.5*(self.Bv / pi*x*self.d_valve)**2) + return Av + + def derivs(self, Core): + """ + Return the position and velocity as an arraym for the valve + + Parameters + ---------- + Core : :class:`PDSimCore ` instance + + Returns + ------- + out_array : :class:`arraym ` instance + + """ + f = empty_arraym(2) + out_array = empty_arraym(2) + x = self.xv[0] + xdot = self.xv[1] + self.get_States(Core) + rho = self.State_up.get_rho() + p_high = self.State_up.get_p() + p_low = self.State_down.get_p() + deltap = (p_high - p_low)*1000 + self.m_vax = self.m_valve*(1 - 0.5*x/self.x_stopper) #Mass of a valve moving as the free body + + if self.valve_type == '1stage': + self.k_valve = 0.11584 *np.exp(x/0.00042) + 762.18456 + + elif self.valve_type == '2stage': + self.k_valve = 0.11584*np.exp(x/0.00042) + 762.18456 + + self.F_fac = 0.04* (x/self.d_port-1)**2 + 0.06 #proportional factor based on the experiment + self.omega_n = np.sqrt(self.k_valve/self.m_vax) #Natural frequency + Dl = self.F_fac*self.Bv*deltap/self.m_vax #Acceleration of the gas againt the valve + + self._pressure_dominant(f,x,xdot,Dl,deltap) + + omega = Core.omega + + out_array[0] = f[0]/omega + out_array[1] = f[1]/omega + + if abs(x) < 1e-15 and xdot < -1e-12: + # print 'stationary valve' + out_array[0] = 0.0 + out_array[1] = 0.0 + + return out_array #[dxdtheta, d(xdot)_dtheta] + + def __dict__(self): + + # print 'Im in dict!' + items=['A_port','A_valve','d_valve','d_port','m_valve','C_D', + 'k_valve','x_stopper','key_up','key_down','xv','State_up', 'State_down'] + + return {item:getattr(self,item) for item in items} + + def __repr__(self): + """ + Return a representation of Valve Model for outputting to screen + """ + # print 'Im in repr!' + s='' + for item in self.__dict__(): + s += item+' : '+str(getattr(self,item))+'\n' + return s + + def __reduce__(self): + # print 'Im in reduce!' + return rebuildValveModel,(self.__dict__(),) + + +######################################################## +# Add isentropicNozzle model into valve model +######################################################## + +def IsentropicNozzle(A, State_up, State_down, other_output = -1): + """ + The mass flow rate is calculated by using isentropic flow model + + Parameters + ---------- + + A : double + Throat area of the nozzle [m\ :math:`^2`\ ] + State_up : :class:`State ` instance + Upstream ``State`` instance + State_down : :class:`State ` instance + Downstream ``State`` instance + other_output : int + Default is to return the mass flow rate, can over-ride by passing ``flow_models.OUTPUT_VELOCITY`` or ``flow_models.OUTPUT_MA`` instead + + Returns + ------- + out : double + Default is to return the mass flow rate, can over-ride by passing flags in the other_output variable + + """ + + #some cython type declarations + # cython.declare(T_up = cython.double, T_down = cython.double, mdot = cython.double) + # Since ideal, R=cp-cv, and k=cp/cv + cp = State_up.get_cp0() + R = 8314.472/State_up.get_MM() #[J/kg/K] + cv = cp-R/1000.0 #[kJ/kg/K] + k = cp / cv + + p_up = State_up.get_p() + T_up = State_up.get_T() + p_down = State_down.get_p() + + # Speed of sound + c = (k*R*T_up)**0.5 + # Upstream density + rho_up = p_up*1000.0/(R*T_up) + pr = p_down/p_up + pr_crit = (1+(k-1)/2)**(k/(1-k)) + + if pr > pr_crit: + # Mass flow rate if not choked [kg/s] + mdot = A*p_up*1000.0/(R*T_up)**0.5*(2*k/(k-1.0)*pr**(2.0/k)*(1-pr**((k-1.0)/k)))**0.5 + # Throat temperature [K] + T_down = T_up*(p_down/p_up)**((k-1.0)/k) + # Throat density [kg/m3] + rho_down = p_down*1000.0/(R*T_down) + # Velocity at throat + v = mdot/(rho_down*A) + # Mach number + Ma = v/c + else: + # Mass flow rate if choked + mdot = A*rho_up*(k*R*T_up)**0.5*(1.+(k-1.)/2.)**((1+k)/(2*(1-k))) + # Velocity at throat + v = c + # Mach Number + Ma = 1.0 + + if other_output < 0: + return mdot + # elif other_output == OUTPUT_VELOCITY: + elif other_output == 2: + return v + elif other_output == OUTPUT_MA: + return Ma + + + + \ No newline at end of file diff --git a/PDSim/heattransfer/HTC.py b/PDSim/heattransfer/HTC.py new file mode 100644 index 0000000..ddbf639 --- /dev/null +++ b/PDSim/heattransfer/HTC.py @@ -0,0 +1,138 @@ +from __future__ import division, print_function +from math import pi, cos, sin, log +import os, sys + +#Math Package +from scipy import optimize +from math import ceil +from scipy.special import erf +from PDSim.misc.scipylike import trapz +from scipy.integrate import quad,simps +from scipy.interpolate import interp1d,splrep,splev +import pylab +import numpy as np +# Coolprop +from CoolProp import State,AbstractState +from CoolProp import State +from CoolProp import CoolProp as CP + + +def HTC(HTCat,T_wall,T_inf,P_film,Ref,L,D_pipe=None,PlateNum=None): + + """ + Nusselt number for different heat transfer categories; + find heat transfer coefficient based on Nusselt number; + characteristic length, A/P; + + Return heat transfer rate by radiation + + Parameters + ---------- + A [m^2] : heat transfer area + epsilon : surface emissivity + T_wall [K]: surface temperature + P_film [kPa]: pressure at the film + Tinf [K]: surrounding temperature + sigma [W/(m^2 K)]: Stefan-Boltzmann constant 5.670*10e-8 + ---------- + Return + h [kW/m^2 K]: Heat transfer coefficient + """ + g = 9.81 # gravity acceleration [m/s^2] + + # Film temperature, used to calculate thermal propertiers + T_film = (T_wall + T_inf)/2 # [K] + + # thermal expansion coefficient, assumed ideal gas + beta = 1/T_film # [1/K] + + # Transport properties calcualtion film + StateFilm = State.State(Ref,dict(T=T_film, P=P_film)) # use the film temperature to find the outer convective coefficient + + Pr_film = StateFilm.Prandtl #[-] + rho_film = StateFilm.rho #[kg/m3] + k_film = StateFilm.k #[kW/m-K] conductivity + mu_film = StateFilm.visc #[Pa-s] + nu_film = mu_film/rho_film # [m^2/s] + cp_film = StateFilm.cp # [kJ/kg/K] + + if T_wall 1e9: + Nu = (0.825 + 0.387*RaL**(1/6) / ( 1 + (0.492/Pr_film)**(9/16) )**(8/27))**2 + else: + Nu = 0.68 + 0.670*RaL**(1/4) / ( 1 + (0.492/Pr_film)**(9/16) )**(4/9) + + elif HTCat == 'horizontal_plate': + if PlateNum == 'upper_surface': + if T_wall > T_inf: # hot plate + if RaL >= 1e4 and RaL <= 1e7: + Nu = 0.54*RaL**(1/4) + elif RaL >= 1e7 and RaL <= 1e11: + Nu = 0.15*RaL**(1/3) + else: + Nu = 0.71*RaL**(1/4) # not quite sure about this equation + else: # cold plate + if RaL >= 1e5 and RaL <= 1e10: + Nu = 0.27*RaL**(1/4) + else: + Nu = 0.71*RaL**(1/4) # not quite sure about this equation, M. AL-ARABI and M. K. EL-RIEDYt; + + elif PlateNum == 'lower_surface': + if T_wall > T_inf: # hot plate + if RaL >= 1e5 and RaL <= 1e10: + Nu = 0.27*RaL**(1/4) + else: + Nu = 0.25*RaL**(1/4) # not quite sure about this equation, M. AL-ARABI and M. K. EL-RIEDYt; + + else: # cold plate + if RaL >= 1e4 and RaL <= 1e7: + Nu = 0.54*RaL**(1/4) + elif RaL >= 1e7 and RaL <= 1e11: + Nu = 0.15*RaL**(1/3) + else: + Nu = 0.25*RaL**(1/4) # not quite sure about this equation, M. AL-ARABI and M. K. EL-RIEDYt; + else: + raise Exception('PlateNum must be either upper_surface or lower_surface') + + elif HTCat == 'horizontal_pipe': + if RaD <= 1e12: + # Churchill and Chu, 1975, RaL->RaD + Nu = (0.60+0.387*RaD**(1/6)/(1 + (0.559/Pr_film)**(9/16))**(8/27))**2 + + else: # Kuehn and Goldstein, 1976. + temp= (( 0.518*(RaD**0.25)*(1+(0.559/Pr_film)**0.6)**(-5/12) )**15 + (0.1*RaD**(1/3))**15)**(1/15) + Nu = 2/(log(1 + 2/temp)) + + elif HTCat == 'vertical_pipe': + if (D_pipe/L) < 35/Gr**(1/4): + F = 1/3*((L/D_pipe)/(1/Gr))**(1/4)+1 + else: + F = 1.0 + Nu = F*(0.825 + 0.387*RaL**(1/6)/(1+(0.492/Pr_film)**(9/16))**(8/27))**2 + + else: + raise Exception('not implemented') + + # convective heat transfer coefficient + if HTCat == 'horizontal_pipe': + h = Nu*k_film/D_pipe + else: + h = Nu*k_film/L + + return h \ No newline at end of file diff --git a/PDSim/heattransfer/__init__.py b/PDSim/heattransfer/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/PDSim/rolling/__init__.py b/PDSim/rolling/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/PDSim/rolling/core.py b/PDSim/rolling/core.py new file mode 100644 index 0000000..986b92a --- /dev/null +++ b/PDSim/rolling/core.py @@ -0,0 +1,1307 @@ +from __future__ import division, print_function, absolute_import +from math import pi, cos, sin, asin,tan,sqrt,log,log10 +from scipy.constants import g +from scipy.integrate import quad,simps +import numpy as np +import os, sys + +from PDSim.flow.flow import FlowPath +from PDSim.flow import flow_models +from PDSim.misc.datatypes import arraym,listm +from PDSim.core.containers import ControlVolume, Tube +from PDSim.core.core import PDSimCore +from PDSim.heattransfer.HTC import HTC +from PDSim.plot.plots import debug_plots +from PDSim.flow.flow_models import ValveModel +from PDSim.rolling import rolling_geo + +#Imports from CoolProp (fluid property database) +from CoolProp import State,AbstractState +from CoolProp import CoolProp as CP +import CoolProp + +class struct(object): + pass + +class RollingPiston(PDSimCore): + + def __init__(self): + #Initialize the base class that RollingPiston is derived from + PDSimCore.__init__(self) + + ## Define the geometry structure + self.geo=rolling_geo.geoVals() + + def set_rolling_geo(self,e,Rc,Rc_outer,Rr,Rs,rv,Hc,b,hv,offsetRoller,D_sp,D_dpc,theta_sm,theta_dm,delta_theta_s,delta_theta_d): + + self.geo.e = e + self.geo.Rc = Rc + self.geo.Rc_outer = Rc_outer + self.geo.Rr = Rr + self.geo.Rs = Rs + self.geo.rv = rv + self.geo.Hc = Hc + self.geo.b = b + self.geo.hv = hv + self.geo.offsetRoller = offsetRoller + self.geo.theta_sm = theta_sm + self.geo.theta_dm = theta_dm + self.geo.D_sp = D_sp + self.geo.D_dpc = D_dpc + self.geo.theta_sm = theta_sm + self.geo.theta_dm = theta_dm + self.geo.delta_theta_s = delta_theta_s + self.geo.delta_theta_d = delta_theta_d + + #Set the flags to ensure all parameters are fresh + self.__Setrolling_geo__=True + + def V_dV_ParkYC(self,theta): + + """ + Park Y.C. 'Transient analysis of a variable speed rotary compressor',Energy Conversion and Management 2010 + Mathematical description of compression chamber + """ + a = self.Rr/self.Rc + Vclearance = 1e-7 #m^3 + V = Vclearance +self.Hc/2*self.Rc**2*((1-a**2)*theta - (1-a)**2*(sin(2*theta))/2 - a**2*asin( (1/a-1)*sin(theta)) -a*(1-a)*sin(theta)*sqrt(1-(a-1)**2*(sin(theta))**2)) - self.b/2*self.Hc*self.Rc*( 1-(1-a)*cos(theta)- sqrt((1-a)**2*(cos(theta))**2+2*a-1) ) + + dVdtheta = (self.Hc*self.Rc**2/2)* (-((1/a-1)*a**2*cos(theta))/(sqrt(1-(1/a-1)**2*(sin(theta))**2)) - a**2 + (1-a)**2*(-cos(2*theta))-a*(1-a)*cos(theta)*sqrt(1-(1/a-1)**2*(sin(theta))**2)+ ((1/a-1)**2*a*(1-a)*(sin(theta))**2*cos(theta) )/(sqrt(1-(1/a-1)**2*(sin(theta))**2)) + 1 ) - (self.b*self.Hc*self.Rc/2)*( (1-a)*sin(theta)+ ((1-a)**2*sin(theta)*cos(theta))/(sqrt((1-a)**2*(cos(theta))**2+2*a-1)) ) + + if V < 0.: + raise Exception('V negative %f' % (V)) + V = 0.0 + return V,dVdtheta + + def Vs_dVs_ParkYC(self,theta): + #Suction chamber + if theta >= 0.0 and theta <= 2*pi: + V,dVdtheta = self.V_dV_ParkYC(theta) + return abs(V),abs(dVdtheta) + elif theta > 2*pi and theta <= 4*pi: + V,dVdtheta = self.V_dV_ParkYC(4*pi - theta) + return abs(V),-abs(dVdtheta) + + def Vc_dVc_ParkYC(self,theta): + #Compression chamber + if theta >= 0.0 and theta <= 2*pi: + V,dVdtheta = self.V_dV_ParkYC(2*pi - theta) + return abs(V),-abs(dVdtheta) + elif theta > 2*pi and theta <= 4*pi: + V,dVdtheta = self.V_dV_ParkYC(theta-2*pi) + return abs(V),abs(dVdtheta) + + def Vs_dVs_ParkYC_2stage(self,theta): + #Second cylinder suction chamber + if theta >= 0.0 and theta <= pi: + V,dVdtheta = self.V_dV_ParkYC(theta + pi) + return abs(V),abs(dVdtheta) + elif theta > pi: + V,dVdtheta = self.V_dV_ParkYC(theta - pi) + return abs(V),abs(dVdtheta) + + def Vc_dVc_ParkYC_2stage(self,theta): + #Second cylinder compression chamber + if theta >= 0.0 and theta <= pi: + V,dVdtheta = self.V_dV_ParkYC(2*pi - theta - pi) + return abs(V),-abs(dVdtheta) + elif theta > pi: + V,dVdtheta = self.V_dV_ParkYC(2*pi - theta + pi) + return abs(V),-abs(dVdtheta) + + def x_theta(self,x): + """ + Vane displacement with the crankshaft angle + """ + a = self.Rr/self.Rc + return self.Rc*(1-(1-a)*cos(x) - sqrt((1-a)**2*(cos(x))**2+2*a-1 ) ) + + def xv_dot(self,theta, stage=1): + """ + Function for calculating velocity of the vane. + """ + R_r = self.Rr + R_v = self.rv + e = self.e + z = 1 - (e*sin(theta)/(R_r+R_v))**2 + alphadot = e*cos(theta)*self.omega/(sqrt(z)*(R_r+R_v)) + + return e*sin(theta)*(alphadot+self.omega) + + def V_disp(self,alpha_s,alpha_d): + """ + Zheng et al. 'Experimental verification of a rolling-piston expander that + applied for low-temperature ORC', Applied Energy 112(2013)1265-1274. + + ParkY.C. 'Transient analysis of a variable speed rotary compressor' + Energy Conversion and Management 51(2010) 277-287 + + """ + V_theta_s = self.Vs_dVs_ParkYC(alpha_s)[0] + V_theta_d = self.Vc_dVc_ParkYC(2*pi - alpha_d)[0] + Int_x_dtheta, err = quad(self.x_theta,0,2*pi) + Vdisp = pi*self.Hc*(self.Rc**2 - self.Rr**2) - self.b*self.Hc*Int_x_dtheta - (V_theta_s + V_theta_d) + + if hasattr(self.mech,'cylinder_config') and self.mech.cylinder_config == 'dual-cylinder': + return 2*Vdisp + elif hasattr(self.mech,'cylinder_config') and self.mech.cylinder_config == 'single-cylinder': + return Vdisp + else: + return Vdisp + + def V_shell(self, theta): + #Approximation of shell volume + V_shell_tot = pi*self.Rshell**2*self.Hshell + Vcylinder = pi*self.Rc_outer**2*self.Hc + Vmotor = pi*self.Rshell**2*self.Hc/3 + V = (V_shell_tot - Vcylinder - Vmotor) + dV = 0.0 + return V,dV + + def xv_simplified(self,theta, stage = 1): + """ + Distance that the vane extends into the cylinder + Assumption: rv = 0.0 + """ + if stage == 2: + R_c = self.Rc_2 + R_r = self.Rr_2 + e = self.e_2 + elif stage == 1: + R_c = self.Rc + R_r = self.Rr + e = self.e + + return R_c - R_r * cos(asin(e / R_r * sin(theta))) - E * cos(theta) + + def A_vb(self,theta, geo): + """ + Leakage area through vertical clearance between vane and cylinder + """ + return 2*self.delta_vb*self.x_theta(theta,geo) + + def A_vt(self): + """ + Leakage area through radial clearance between vane tip and roller + """ + return self.hv*self.delta_vt + + def A_f(self): + """ + Leakage area between compression chamber and muffler + """ + return pi/4*self.D_muff**2 + + def A_lc(self): + """ + Leakage area between upper and lower shell + """ + return self.A_ls + + def A_uc(self): + """ + Leakage area from discharge port + """ + return pi/4*self.D_uc**2 + + def A_rc(self,theta, delta): + """ + Function for calculating area of leakage path across roller. + """ + x = self.e*sin(theta)/self.Rr + alpha = asin(x) + beta = alpha + theta + + return 2*(2*pi - beta)*self.Rr_i*delta + + def A_10(self): + """ + Leakage area between accumulator and suction pipe + """ + return self.A_sp + + def A_21(self,theta, stage=1): + """ + Leakage area between suction chamber and suction pipe + """ + l_theta_max = pi*self.D_sp + + if theta <= self.theta_s1: + l_theta = 0.0 + A21 = 0.0 + elif theta > self.theta_s1 and theta <= self.theta_sm: + l_theta = l_theta_max/2*(theta - self.theta_s1)/(self.theta_sm - self.theta_s1) + phi = theta - self.theta_s1 + z_theta = self.xv_simplified(phi, stage) + A21 = z_theta*l_theta/2 + elif theta > self.theta_sm and theta <= self.theta_s2: + l_theta = l_theta_max*(1 - 0.5*(self.theta_s2 - theta)/(self.theta_sm - self.theta_s1)) + phi = theta - self.theta_s1 + z_theta = self.xv_simplified(phi, stage) + A21 = z_theta*l_theta/2 + elif theta > self.theta_s2: + l_theta = l_theta_max + phi1 = theta - self.theta_s1 + phi2 = theta - self.theta_s2 + z_theta = (self.xv_simplified(phi1, stage) + self.xv_simplified(phi2, stage))/2 + A21 = z_theta*l_theta + + if A21 > self.A_sp: + A21 = self.A_sp + + return A21 + + def A_31(self,theta, stage=1): + """ + Leakage area between compression chamber and suction pipe + """ + l_theta_max = pi*self.D_sp + + if theta <= self.theta_s1: + l_theta = l_theta_max + elif theta > self.theta_s1 and theta <= self.theta_sm: + l_theta = l_theta_max*(1 - 0.5*(theta-self.theta_s1)/(self.theta_sm-self.theta_s1)) + elif theta > self.theta_sm and theta < self.theta_s2: + l_theta = l_theta_max/2*(self.theta_s2 - theta)/(self.theta_sm - self.theta_s1) + elif theta >= self.theta_s2: + l_theta = 0 + + phi = self.theta_s2 - theta + z_theta = self.xv_simplified(phi, stage) + A31 = z_theta*l_theta/2 + + if A31 > self.A_sp: + A31 = self.A_sp + + return A31 + + def A_32(self,theta, stage=1): + """ + Flank leakage area between compression chamber and suction chamber + """ + if theta <= self.theta_d1: + H32 = self.Hc + elif theta > self.theta_d1 and theta < self.theta_dm: + H32 = self.Hc - self.hp*(theta-self.theta_d1)/(self.theta_dm-self.theta_d1) + else: + H32 = self.Hc - self.hp*(self.theta_d2-theta)/(self.theta_dm-self.theta_d1) + + if H32 > self.Hc: + H32 = self.Hc + + return H32*self.delta_gap + + def A_42_ref(self,stage=1): + """ + Reference value used to calculate flow coefficient from clearance volume to suction chamber + """ + + phi = self.theta_d2 - self.theta_d1 + z_theta = self.xv_simplified(phi, stage) + + return self.l_dm*z_theta/2 + + def A_42(self,theta, stage=1): + """ + Leakage area between clearance volume and suction chamber + """ + phi = self.theta_d2 - self.theta_d1 + z_theta = self.xv_simplified(theta, stage) + + if theta < self.theta_d1: + A42_side = 0.0 + A42_t4 = 0.0 + A42 = 0.0 + elif theta >= self.theta_d1 and theta < self.theta_d2: + l_theta = self.l_dm * (self.theta - self.theta_d1)/(self.theta_d2 - self.theta_d1) + phi = theta - self.theta_d1 + z_theta = self.xv_simplified(phi, stage) + A42_side = l_theta*z_theta/2 + A42_t4 = z_theta*(self.D_dp- self.D_dpc)/2 + A42 = A42_side + A42_t4 + elif theta >= self.theta_d2: + phi1 = self.theta_d2 - self.theta_d1 + phi2 = theta - self.theta_d1 + z_theta1 = self.xv_simplified(phi1, stage) + z_theta2 = self.xv_simplified(phi2, stage) + A42_t4 = z_theta2*(self.D_dp - self.D_dpc)/2 + A_42 = self.l_dm*z_theta1/2 # + A42_t4 + + return A42 + + def A_43_ref(self,theta, stage=1): + """ + Reference value used to calculate flow coefficient from clearance volume to compression chamber + """ + A43a_max = pi/8*self.D_dp**2 + phi1 = self.theta_d1 - theta + phi2 = self.theta_d2 - theta + phim = self.theta_dm - theta + + z_thetam = self.xv_simplified(phim, stage) + if z_thetam > self.D_dp/2: + z_thetam = self.D_dp/2 + + if theta < self.theta_d1: + z_theta2 = (self.xv_simplified(phi1, stage) + self.xv_simplified(phi2, stage))/2 + else: + z_theta2 = self.xv_simplified(phi2, stage)/2 + + A43_side = z_thetam*self.hp/4 + A43_a = self.D_dp*z_theta2 + if A43_a > A43a_max: + A43_a = A43a_max + A43_b = pi*self.D_dp*self.hp/4 + A43_ref = sqrt(A43_a**2 + A43_b**2) + A43_side + + return A43_ref + + def A_43(self,theta, stage=1): + """ + Leakage area between clearance volume and compression chamber + """ + A43a_max = pi/8*self.D_dp**2 + phi1 = self.theta_d1 - theta + phi2 = self.theta_d2 - theta + phim = self.theta_dm - theta + + z_thetam = self.xv_simplified(phim, stage) + if z_thetam > self.D_dp/2: + z_thetam = self.D_dp/2 + + if theta < self.theta_d1: + z_theta2 = (self.xv_simplified(phi1, stage) + self.xv_simplified(phi2, stage))/2 + else: + z_theta2 = self.xv_simplified(phi2, stage)/2 + + l_theta = self.l_dm*(self.theta_d2 - theta)/(self.theta_d2 - self.theta_d1) + if l_theta < 0.0: + l_theta = 0.0 + elif l_theta > self.l_dm: + l_theta = self.l_dm + + A43_side = z_thetam*self.hp/4 + A43_a = self.D_dp*z_theta2 + if A43_a > A43a_max: + A43_a = A43a_max + + A43_b = pi*self.D_dp*self.hp/4 + A43_c1 = sqrt(A43_a**2 + A43_b**2) + A43_side + A43_c2 = l_theta*z_theta2 + + if A43_c1 > A43_c2: + A43 = A43_c2 + else: + A43 = A43_c1 + + return A43 + + def calculate_heat_transfer_area(self,theta): + """ + Calculation of chamber total heat transfer area + """ + a = self.Rr/self.Rc + hv_theta = self.Rc*(1-(1-a)*cos(theta)-sqrt((1-a)**2*(cos(theta))**2) +2*a-1) + f_theta = (1-a**2)*theta-0.5*(1-a)**2*sin(2*theta) - a**2*asin((1/a-1)*sin(theta)) - a*(1-a)*sin(theta)*sqrt(1-(1/a-1)**2*(sin(theta))**2) + phi = pi/2 + theta - asin((self.Rc - hv_theta + self.e*cos(pi - theta))/self.Rr ) + + #delta = self.Rc- (2*self.e*cos(theta-phi)+sqrt(4*e**2*(cos(theta-phi))**2-4*(e**2 - self.Rr**2)))/2 + + A = theta*(self.Rc + self.Rr)*self.Hc+ 2*0.5*self.Rc**2*f_theta + + if A < 0.0: + raise Exception('A_ht < 0.0') + + return A + + def A_HT_s(self,theta): + #Suction chamber + if theta >= 0.0 and theta <= 2*pi: + A = self.calculate_heat_transfer_area(theta) + return A + elif theta > 2*pi and theta <= 4*pi: + A = self.calculate_heat_transfer_area(4*pi-theta) + return A + + def A_HT_c(self,theta): + #Compression chamber + if theta >= 0.0 and theta <= 2*pi: + A = self.calculate_heat_transfer_area(2*pi - theta) + return A + elif theta > 2*pi and theta <= 4*pi: + A = self.calculate_heat_transfer_area(theta-2*pi) + return A + + def A_HT_s_2stage(self,theta): + #Second cylinder suction chamber + if theta >= 0.0 and theta <= pi: + A = self.calculate_heat_transfer_area(theta + pi) + return A + elif theta > pi: + A = self.calculate_heat_transfer_area(theta - pi) + return A + + def A_HT_c_2stage(self,theta): + #Second cylinder compression chamber + if theta >= 0.0 and theta <= pi: + A = self.calculate_heat_transfer_area(2*pi - theta - pi) + return A + elif theta > pi: + A = self.calculate_heat_transfer_area(2*pi - theta + pi) + return A + + def r_chamber_mean(self,theta): + """ + Calculates the mean chamber radius to be used for heattransfer coefficient calculation + """ + a = self.Rr/self.Rc + N = 100 + phi_array = np.linspace(0,theta,N) + R_array = np.zeros(len(phi_array)) + ravg_array = np.zeros(len(phi_array)) + + for i in range(len(phi_array)): + R_array[i] = self.Rc*((1-a)*cos(theta-phi_array[i])+sqrt((1-a)**2*(cos(theta- phi_array[i]))**2+(2*a-1))) + ravg_array[i] = (self.Rc + R_array[i])/2 + + r = np.sum(ravg_array)/N + + return r + + def alp(self,p,T): + """ + Solubility Refrigerant in oil + p: pressure [kPa] + T: temperature [K] + """ + if self.Ref == 'R410a' or self.Ref == 'HEOS::R410a' or self.Ref == 'REFPROP::R401a.mix': + #Solubility of the R410A and PVE Oil + P0 = p*1e-3 #pressure [MPa] + T0 = T - 273.15 #Temperaure [C] + + Sol20 = 3.49058 + 14.33071*P0 + 3.94294*P0**2 + Sol40 = 5.57178 + 0.76551*P0 + 5.97508*P0**2 + Sol60 = 4.81862 - 0.34116*P0 + 3.51925*P0**2 + Sol80 = 2.88481 + 2.19068*P0 + 1.54295*P0**2 + Sol100 = 2.2739 + 2.24531*P0 + 1.07309*P0**2 + Sol120 = 0.56494 + 3.15252*P0 + 0.65613*P0**2 + + if T0 < 20: + alp = Sol20 + elif 20 <= T0 and T0 < 40: + alp = Sol40 - (Sol40 - Sol20)*(40 - T0)/(40 - 20) + elif T0 >= 40 and T0 < 60: + alp = Sol60 - (Sol60 - Sol40)*(60 - T0)/(60 - 40) + elif T0 >= 60 and T0 < 80: + alp = Sol80 - (Sol80 - Sol60)*(80 - T0)/(80 - 60) + elif T0 >= 80 and T0 < 100: + alp = Sol100 - (Sol100 - Sol80)*(100 - T0)/(100 - 80) + elif T0 >= 100 and T0 <= 120: + alp = Sol120 - (Sol120 - Sol100)*(120 - T0)/(120 - 100) + elif T0 > 120: + alp = Sol120 + + if alp <= 0: + alp = 0 + if alp > 10: + alp = 10 + + else: + alp = 0.0 + + return alp/100 + + def Oil_viscosity(self,p,T): + """ + Used to calculate oil viscosity as a function of temperature and pressure. + p : pressure kPa + T : temperature K + """ + alp1 = self.alp(p, T)*100 + T1 = T - 273.15 + + if T1 < 60: + #0<=T< 60 + mu_10 = 196.05-13.05438*T1 + 0.39545*T1**2 - 0.00617*T1**3 + 0.0000478516*T1**4 - 0.000000145573*T1**5 + elif T1 >= 60: + #60<=T<=100 + mu_10 = 59.9 - 1.52308*T1 + 0.01471*T1**2 - 0.0000504167*T1**3 + + if T1 < 20: + #-20<=T<20 + mu_20 = 44.75 - 2.83417*T1 + 0.10098*T1**2 - 0.00198*T1**3 + 0.0000212222*T1**4 - 0.000000116146*T1**5 + 0.000000000256076*T1**6 + elif T1 >= 20: + #20<=T<=100 + mu_20 = 24.682 - 0.5437*T1 + 0.00472*T1**2 - 0.0000145833*T1**3 + + if T1 < 60: + #-20<=T<60 + mu_30 = 13.78 - 0.69392*T1 + 0.03069*T1**2 - 0.000867812*T1**3 + 0.0000135104*T1**4 - 0.000000105937*T1**5 + 0.000000000326823*T1**6 + elif T1 >= 60: + #60<=T<=100 + mu_30 = 10.15 - 0.184*T1 + 0.00128*T1**2 - 0.0000025*T1**3 + + mu_40 = 6.60955 - 0.20951*T1 + 0.00385*T1**2 - 0.0000303243*T1**3 + 0.0000000262784*T1**4 + 0.000000000507812*T1**5 + mu_50 = 3.18206 - 0.07856*T1 + 0.00203*T1**2 - 0.0000350712*T1**3 + 0.000000323208*T1**4 - 0.00000000115685*T1**5 + + if alp1 < 20: + mu_oil = (mu_20 - ((20 - alp1)*(mu_20 - mu_10)/10))*0.001 + elif alp1 >= 20 and alp1 < 30: + mu_oil = (mu_30 - ((30 - alp1)*(mu_30 - mu_20)/10))*0.001 + elif alp1 >= 30 and alp1 < 40: + mu_oil = (mu_40 - ((40 - alp1)*(mu_40 - mu_30)/10))*0.001 + elif alp1 >= 40: + #And alp1 <= 90 Then + mu_oil = (mu_50 - ((50 - alp1)*(mu_50 - mu_40)/10))*0.001 + else: + print(alp1, T1) + raise Exception('Out of range') + + return mu_oil + + def Inlet(self, FlowPath, **kwargs): + #if self.theta >= self.psi: + FlowPath.A=self.A_suction + mdot=flow_models.IsentropicNozzle(FlowPath.A,FlowPath.State_up,FlowPath.State_down) + return mdot + + def Outlet(self, FlowPath, **kwargs): + FlowPath.A=self.A_discharge + mdot=flow_models.IsentropicNozzle(FlowPath.A,FlowPath.State_up,FlowPath.State_down) + return mdot + + def Suction(self,FlowPath,**kwargs): + if FlowPath.key_up=='A': + # pressure in compressor higher than the inlet line + # valve is closed - no flow + return 0.0 + if self.theta >= self.theta_s1: + FlowPath.A=self.A_suction + mdot=flow_models.IsentropicNozzle(FlowPath.A,FlowPath.State_up,FlowPath.State_down) + return mdot + else: + return 0.0 + + def Suction_2stage(self,FlowPath,**kwargs): + if FlowPath.key_up=='C': + # pressure in compressor higher than the inlet line + # valve is closed - no flow + return 0.0 + if self.theta <= pi or self.theta >= self.theta_s1+pi: + FlowPath.A=self.A_suction + mdot=flow_models.IsentropicNozzle(FlowPath.A,FlowPath.State_up,FlowPath.State_down) + return mdot + else: + return 0.0 + + def Discharge(self,FlowPath,**kwargs): + if FlowPath.key_down=='B': + # pressure in compressor lower than the discharge line + # valve is closed - no flow + return 0.0 + else: + if self.theta <= self.theta_d2: + FlowPath.A=self.discharge_valve.A() + mdot=flow_models.IsentropicNozzle(FlowPath.A,FlowPath.State_up,FlowPath.State_down) + return mdot + else: + return 0.0 + + def Discharge_2stage(self,FlowPath,**kwargs): + if FlowPath.key_down=='D': + # pressure in compressor lower than the discharge line + # valve is closed - no flow + return 0.0 + else: + if self.theta <= self.theta_d2 - pi: + FlowPath.A=self.discharge_valve_upper.A() + mdot=flow_models.IsentropicNozzle(FlowPath.A,FlowPath.State_up,FlowPath.State_down) + return mdot + else: + return 0.0 + + def CompressibleGasLeakage_vb(self, FlowPath, **kwargs): + """ + From compression to suction chamber through vertical clearance between vane and cylinder + """ + Cv = self.Cflow_vb + FlowPath.A=self.delta_vb*self.b + mdot=Cv*flow_models.IsentropicNozzle(FlowPath.A,FlowPath.State_up,FlowPath.State_down) + return mdot + + def CompressibleGasLeakage_vt(self, FlowPath, **kwargs): + """ + From compression to suction chamber through radial clearance between vane tip and rolling piston + """ + Cv = self.Cflow_vt + FlowPath.A=self.delta_vt*self.Hc + mdot=Cv*flow_models.IsentropicNozzle(FlowPath.A,FlowPath.State_up,FlowPath.State_down) + return mdot + + def CompressibleGasLeakage_32(self, FlowPath, **kwargs): + """ + From compression to suction chamber through radial clearance between cylinder wall and + rolling piston + """ + Cv = self.Cflow_32 + FlowPath.A=self.A_32(self.theta) + mdot=Cv*flow_models.IsentropicNozzle(FlowPath.A,FlowPath.State_up,FlowPath.State_down) + return mdot + + def CouettePoiseuilleLeakage_VaneSlot(self, FlowPath, **kwargs): + State_up = FlowPath.State_up + State_down = FlowPath.State_down + + T_up = State_up.get_T() + P_up = State_up.get_p() + rho_up = State_up.get_rho() + + T_down = State_down.get_T() + P_down = State_down.get_p() + rho_down = State_up.get_rho() + + if FlowPath.key_up=='shell': + # Vane slot chamber + x_up = self.x_slot + if FlowPath.key_down=='A': + # Suction chamber + x_down = self.xs + elif FlowPath.key_down=='B': + # Compression chamber + x_down = self.xc + elif FlowPath.key_up=='B': + # Compression chamber + x_up = self.xc + # Vane slot chamber + x_down = self.x_slot + elif FlowPath.key_up=='A': + # Suction chamber + x_up = self.xs + # Vane slot chamber + x_down = self.x_slot + else: + raise Exception('VaneSlotLeakage condition not defined') + + #Volume flow rate calculated using mixed Plane Couette and Poiseuille flow. + mu_gas_up = State_up.get_visc() + mu_oil_up = self.Oil_viscosity(P_up, T_up) + mu_up = x_up*mu_oil_up + (1 - x_up)*mu_gas_up + mu_gas_down = State_down.get_visc() + mu_oil_down = self.Oil_viscosity(P_down, T_down) + mu_down = x_down * mu_oil_down + (1 - x_down) * mu_gas_down + mu_mean = (mu_up + mu_down)/2 + + xv_dot = self.xv_dot(self.theta) + + Vdot_eb = self.hv*self.delta_vb*(xv_dot/2 + (P_up - P_down)*self.delta_vb**2/(12*mu_mean*self.L_slot)) + + if Vdot_eb >= 0: + rho_mix = 1/(x_up/self.rho_oil + (1 - x_up)/rho_up) + else: + rho_mix = 1/(x_down/self.rho_oil + (1 - x_down)/rho_down) + + mdot = Vdot_eb*rho_mix + + return mdot + + def LaminarViscousLeakage_Roller(self, FlowPath, **kwargs): + """ + Used to calculate mass flow rate for leakage across the roller, modeled as viscous flow. + """ + # Rotation angle + theta = self.theta + + State_up = FlowPath.State_up + State_down = FlowPath.State_down + + T_up = State_up.get_T() #[K] + P_up = State_up.get_p() #[kPa] + rho_up = State_up.get_rho() #[kg/m3] + + T_down = State_down.get_T() #[K] + P_down = State_down.get_p() #[kPa] + rho_down = State_up.get_rho() #[kg/m3] + + if FlowPath.key_up=='shell': + # Oil concentration across roller + x_up = self.x_roller + if FlowPath.key_down=='A': + # Suction chamber + x_down = self.xs + elif FlowPath.key_down=='B': + # Compression chamber + x_down = self.xc + elif FlowPath.key_up=='B': + # Compression chamber + x_up = self.xc + # Oil concentration across roller + x_down = self.x_roller + elif FlowPath.key_up=='A': + # Suction chamber + x_up = self.xs + # Oil concentration across roller + x_down = self.x_roller + else: + raise Exception('RollerLeakage condition not defined') + + mu_gas_up = State_up.get_visc() #[Pa-s] + mu_oil_up = self.Oil_viscosity(P_up, T_up) + mu_up = x_up*mu_oil_up + (1 - x_up)*mu_gas_up + mu_gas_down = State_down.get_visc() + mu_oil_down = self.Oil_viscosity(P_down, T_down) + mu_down = x_down * mu_oil_down + (1 - x_down) * mu_gas_down + mu_mean = (mu_up + mu_down)/2 + + x = self.e*sin(theta)/self.Rr + alpha = asin(x) + beta = alpha + theta + + if FlowPath.key_down=='B' or FlowPath.key_up=='B' : + A = 2*(2*pi - beta)*(self.Rr + self.Rr_i)/2*self.delta_side + elif FlowPath.key_up=='A' or FlowPath.key_down=='A': + A = 2*beta*(self.Rr + self.Rr_i)/2*self.delta_side + + if FlowPath.key_up=='shell' and FlowPath.key_down=='B' : + DeltaP = (P_up - P_down) # [kPa] + Vdot_rc = 2*pi*self.delta_side**3*DeltaP/(6*mu_mean*log(self.Rr/self.Rr_i))*((2*pi - beta)/(2 * pi)) + rho_mix = 1/(x_up/self.rho_oil + (1 - x_up)/rho_up) + mdot = rho_mix*Vdot_rc + else:#FlowPath.key_up=='B' or FlowPath.key_up=='A': + Cv = self.Cflow_roller + mdot=Cv*flow_models.IsentropicNozzle(A,FlowPath.State_up,FlowPath.State_down) + # else: + # raise Exception('RollerLeakage condition not defined') + return mdot + + def TubeCode(self, Tube): + """ + A thin wrapper of the isothermal wall tube from flow_models.py + """ + Tube.Q = flow_models.IsothermalWallTube(Tube.mdot, + Tube.State1, + Tube.State2, + Tube.fixed, + Tube.L, + Tube.ID, + T_wall = self.Tlumps[0]) + + def heat_transfer_callback(self, theta): + """ + A callback used by PDSimCore.derivs to calculate the heat transfer + to the gas in the working chamber. + We return an arraym instance the same length as the number of CV in existence + + """ + if self.HT_on == True: + Qarray = [] + + for key in self.CVs.exists_keys: + + if key == 'A': + #Suction chamber + V,dVdtheta = self.Vs_dVs_ParkYC(theta) + A_ht = self.A_HT_s(theta) #[m2] + + elif key == 'B': + #Compression chamber + V,dVdtheta = self.Vc_dVc_ParkYC(theta) + A_ht = self.A_HT_c(theta) #[m2] + + if key == 'C': + #Suction chamber + V,dVdtheta = self.Vs_dVs_ParkYC_2stage(theta) + A_ht = self.A_HT_s_2stage(theta) #[m2] + + elif key == 'D': + #Compression chamber + V,dVdtheta = self.Vc_dVc_ParkYC_2stage(theta) + A_ht = self.A_HT_c_2stage(theta) #[m2] + + #Assign lumped temperature to wall + T_w = self.Tlumps[0] #[K] + + #Calculate thermophysical properties + cp = self.CVs[key].State.cp # kJ/kg/K + rho = self.CVs[key].State.rho #[kg/m3] + k = self.CVs[key].State.k #[kW/m-K] + mu = self.CVs[key].State.visc #[Pa-s] + T = self.CVs[key].State.T #[K] + + #Avoid division by zero when Aht or Dh are zero + r_avg = (self.Rc + self.Rr)/2 # self.r_chamber_mean(theta) + if A_ht > 0.0 and V > 0.0: + D_h = 4*V/A_ht # Liu (1994) Eq. 3.4.2 [m] + if D_h < 0.0: + raise Exception('D_h < 0.0') + + if D_h > 0.0: + u = D_h*self.omega/2 # Liu (1994) Eq. 3.4.3 #[m/s] + Pr = (cp*mu)/k #[-] + Re = (rho*u*D_h)/mu #[kg/m3*m/s*m/Pa/s]=[-] + + h_c = 0.025*(k/D_h)*Pr**(0.4)*Re**(0.8)*(1.0+1.77*(D_h/r_avg)) #[kW/m2/K] Liu (1994) Eq. 3.4.8 + Q = h_c*A_ht*(T_w-T) #Net heat into control volume [kW] + else: + Q = 0.0 #kW + else: + Q = 0.0 #kW + + if key == 'shell': + T_w = self.Tlumps[0] # [K] + T = self.CVs[key].State.T # [K] + P_film = self.CVs[key].State.p # [kPa] + + # Vertical side shell + h_c_vertical = HTC('vertical_plate',T_w,T,P_film,self.Ref,self.Hshell) + A_shell_vertical = 2*pi*self.Rshell*self.Hshell # m2 + Q_vertical = h_c_vertical*A_shell_vertical*(T_w-T) + + # Top Surface + h_c_topshell = HTC('horizontal_plate',T_w,T,P_film,self.Ref,self.Rshell*2,PlateNum='upper_surface') + A_shell_top = pi*self.Rshell**2 # m2 + Q_topshell = h_c_topshell*A_shell_top*(T_w-T) + + # Bottom Surface + h_c_bottomshell = HTC('horizontal_plate',T_w,T,P_film,self.Ref,self.Rshell*2,PlateNum='lower_surface') + A_shell_bottom = pi*self.Rshell**2 # m2 + Q_bottomshell = h_c_bottomshell*A_shell_bottom*(T_w-T) + + # Total heat transfer + Q = Q_vertical + Q_topshell + Q_bottomshell + + Qarray.append(Q) + + return arraym(Qarray) + + elif self.HT_on == False: + return arraym([0.0]*self.CVs.N) + + else: + raise Exception('HT_on is either True or False') + + + def calculate_force_terms(self,p_shell = None): + """ + Calculate the force profiles, mean forces, moments, etc. + + Parameters + ---------- + p0 : float, or class instance + If a class instance, must provide a function __call__ that takes as its first input the Rolling class + + """ + if not isinstance(p_shell,float): + raise ValueError('calculate_force_terms must get a float top vane pressure for now') + + Nstep = self.Itheta+1 + _slice = range(self.Itheta+1) + theta = self.t[_slice] + + e = self.e + omega = self.omega + epsilon = e/self.Rr + hv = self.hv #vane height + b = self.b #vane width + rv = self.rv + Rr = self.Rr + Re = self.Re + Rs = self.Rs + Hc = self.Hc + rho_v = 7200 #kg/m3 + rho_r = 7200 #kg/m3 + I_r = pi*rho_r*(Rr**4 - Re**4)*Hc/2 #Moment of inertia of rolling piston + + #Volume of the vane from Liu PhD Thesis + av = rv - sqrt(rv**2 - b**2/4 ) + theta_v = asin(b/(2*rv)) + A5 = rv**2*theta_v- b/2*sqrt(rv**2 - b**2/4) + Av = b*(hv-av) + Vv = Hc*(Av + A5) + mv = rho_v*Vv + + #vane displacement + xv = e*(1-np.cos(theta)) + Rr*(1 - np.sqrt(1-epsilon**2*(np.sin(theta))**2)) + #vane sliding velocity + xdot_v = e*omega*(np.sin(theta) + (epsilon*np.sin(2*theta))/(2*np.sqrt(1-epsilon**2*(np.sin(theta))**2))) + #vane acceleration + self.xdotdot_v = e*omega**2*(np.cos(theta) + (epsilon*np.cos(2*theta))/(np.sqrt(1- epsilon**2*(np.sin(theta))**2)) + (epsilon**3*np.sin(2*theta))/(4*(np.sqrt( (1-epsilon**2*(np.sin(theta))**2)**3 ))) ) + alpha = np.arcsin(e*np.sin(theta)/(Rr+rv)) + pc = self.p[1][_slice]*1000 #Pa + pe = self.p[0][_slice]*1000 #Pa + + Tc = self.T[1][_slice] #K + Te = self.T[0][_slice] #K + + #Forces on the vane + self.Fc = Hc*(b*p_shell - pe*(b/2 - rv*np.sin(alpha)) - pc*(b/2 + rv*np.sin(alpha))) + self.Fh = Hc*xv*(pc-pe) + self.F_lv = -mv*self.xdotdot_v + x0 = 0.0 #initial extension spring + K_v = 1150 #spring stiffness [N/m] can be obtained from Ooi 1994 800-3000 N/m + self.Fk = K_v*(x0 - np.fabs(xv)) + + #Forces on the roller + m_r = rho_r*(pi*(Rr**2-Re**2)*Hc) + #Centrifugal force on roller + self.F_lp = m_r*e*omega**2 + + """ + From Yanagisawa and Shimizu "Friction losses in rolling piston type compressor III" + mu_v = 0.15 - 35*np.sqrt(mu_v/(Fn/Hc)) + """ + mu_v = 0.16 + mu_s = 0.18 #typically 0.16-0.22 + delta_rv = rv*(1-np.cos(alpha)) + self.Fn = (mu_s*self.Fh*(hv+mu_s*b)+(self.Fk+self.Fc+self.F_lv)*(xv-hv))/( mu_s*(np.sin(alpha)-mu_v*np.cos(alpha))*(xv+hv+mu_s*b-2*delta_rv) + (np.cos(alpha) + mu_v*np.sin(alpha))*(xv - hv - 2*mu_s*rv*np.sin(alpha)) ) + + self.F_R1 = 0.5*self.Fh + 1/(2*mu_s)*(self.Fk + self.Fc + self.F_lv) + 0.5*self.Fn*(mu_v*np.cos(alpha) - (mu_v/mu_s)*np.sin(alpha) - np.sin(alpha) - np.cos(alpha)/mu_s) + self.F_R2 = -0.5*self.Fh + 1/(2*mu_s)*(self.Fk + self.Fc + self.F_lv) - 0.5*self.Fn*(mu_v*np.cos(alpha) + (mu_v/mu_s)*np.sin(alpha) - np.sin(alpha) + np.cos(alpha)/mu_s) + + #Friction force generated by radial force of vane on roller + self.Ft = mu_v*self.Fn + self.Ft_R1 = mu_s*self.F_R1 + self.Ft_R2 = mu_s*self.F_R2 + # force caused by pressure difference between suction and discharge chambers + self.Fg = 2*Rr*Hc*(pc-pe)*np.sin((theta+alpha)/2) + #Angular position of Fg + self.theta_Fg = (theta - alpha)/2 + + #Force balance equations on roller + self.Fr = self.Fg*np.cos((theta+alpha)/2) - self.Fn*np.cos(theta+alpha) - self.Ft*np.sin(theta+alpha) + self.F_lp + self.Ftheta = -self.Fg*np.sin((theta+alpha)/2) + self.Fn*np.sin(theta+alpha) - self.Ft*np.cos(theta+alpha) + + #Resultant + self.F = np.sqrt(self.Fr**2 + self.Ftheta**2) + self.theta_f = theta + np.arctan(self.Ftheta/self.Fr) + + #Gas torque + self.Mg = e*self.Fg*np.sin((theta+alpha)/2) + #Torque generated by radial force on the center of roller imposed by vane + self.Mn = -e*self.Fn*np.sin(theta+alpha) + #Friction torque generated by mu_v*Fn = Ft + #TODO: check sign + self.Mn_t = e*self.Ft*np.cos(theta+alpha) + #Moment acting on the crankshaft by the roller + self.Mf = self.Mg + self.Mn + self.Mn_t + + #Assming steady-state conditions + Toil = (Tc+Te)/2 #K + rho_oil = -0.0034*(Toil-273.15)**2-0.0776*(Toil-273.15) + 960.88 + #Vogel equation + mu_oil = rho_oil*(1.742e-4)*np.exp(3664/(Toil-10.27))*1e-6 #Toil [K] + #Clearences + delta_e = 0.25e-3 #[m] + delta_1 = 0.25e-3 #[m] + delta_2 = 0.25e-3 #[m] + delta_3 = 0.25e-3 #[m] + delta_4 = 0.25e-3 #[m] + + #Simplified analysis roller rotational speed + l_r = np.sqrt(Re**2+e**2-2*Re*e*np.cos(pi-self.theta_f+theta)) + cos_thetar = np.sqrt(1 - (e/l_r*np.sin(self.theta_f - theta))**2) + self.omega_r = (2*pi*mu_oil*omega*Hc*Re**3/delta_2 - Rr*self.Ft)*delta_2*delta_e/(2*pi*mu_oil*(Hc*Re**3 + (Hc*Re**3*delta_e)*delta_2)) + +# self.omega_r = (mu_oil*omega*Hc*Re**3/delta_2 - Rr*self.Ft/(2*pi))/(mu_oil*Hc*Re**3/delta_2 + mu_oil*(Rr**4 - Re**4)/delta_1) + + #Relative angular velocity + DELTA_omega = omega - self.omega_r + #Frictional loss between eccentric and the inner surface of the roller Ler + self.Ler = DELTA_omega*(2*pi*mu_oil*DELTA_omega)*Hc*Re**3/delta_2 + + #Frictional loss between the roller face and the cylinder head faces Lrc + self.Lrc = omega*(2*pi*self.omega_r*mu_oil*(Rr**4-Re**4))/delta_1 + + #Frictional loss between eccentric face and the cylinder head faces Lec + self.Lec = omega*(pi*omega*mu_oil*(2*Re**4-Rs**4))/delta_3 + + #Frictional loss between vane tip and roller Lv + self.Lv = self.Ft*(Rr*self.omega_r + (e*omega*np.cos(theta))/np.cos(alpha)) + + #Frictional loss between vane sides and vane slot Ls + alphadot = 1/omega*(e*np.cos(theta)/(Rr+rv))/(np.sqrt(1 - (e*np.sin(theta)/(Rr+rv))**2)) + xdot = e*omega*np.sin(theta) + (Rr+rv)*alphadot*np.sin(alpha) + self.Ls = (self.Ft_R1 + self.Ft_R2)*np.fabs(xdot) + #Frictional loss between shaft and upper journal bearing Lb2 + h_mj2 = 0.067 # length bearing[m] + self.Lbs = omega*(2*pi*mu_oil*omega*h_mj2*Rs**3)/delta_4 + + #Frictional loss between shaft and lower bearing , Lb1 + h_sj = 0.025 # length bearing[m] + self.Lbl = omega*(2*pi*mu_oil*omega*h_sj*Rs**3)/delta_4 + + #Lower and upper balancer distance and eccentricity + l_b1 = 0.1241 #[m] + l_b2 = 0.2221 #[m] + + #Eccentric force distance + l_ec = 0.399 #[m] + + #Mass crank + rho_crank = 8000 #kg/m3 + l_c = 0.200 #[m] + m_ec = rho_crank*(pi/2*(Re-Rs)**2*l_c) + + F_ec = (m_r+m_ec)*e*omega**2 + F_b2 = F_ec*(l_ec-l_b1)/(l_b1-l_b2) + F_b1 = F_ec + F_b2 + + #Estimation of balancers + e_b1 = 0.03 #[m] + m_b1 = e/e_b1*(m_r + m_ec)*(1+ (l_b1-l_ec)/(l_b2-l_b1)) + e_b2 = 0.03 #[m] + m_b2 = e/e_b2*(m_r + m_ec)*(l_b1-l_ec)/(l_b2-l_b1) + + def friction_losses(self,shell_pressure = 'mean:shell'): + + """ + Friction losses are caused by rubbing parts, bearings etc + + Ooi K.T. 'Design optimization of a rolling piston compressor + for refrigeration', Applied Therm. Eng. 25(2005) 813-829 + + Parameters + ---------- + shell_pressure : string, 'low', 'low:shell', 'mid', or 'high' + + low uses the upstream pressure of the machine, + + low:shell uses the pressure after the inlet tube + + mid uses the average of upstream and downstream pressures + + high uses the pressure downstream of the machine + + """ + #inlet pressure [Pa] + inlet_pressure = self.Tubes.Nodes[self.key_inlet].p*1000 + outlet_pressure = self.Tubes.Nodes[self.key_outlet].p*1000 + + # Find the tube with the inlet node + Tube = self.Tubes[self.key_inlet] + # Get the state that is not the inlet state + if Tube.key1 == 'self.key_inlet': + shell_pressure_val = Tube.State2.p + else: + shell_pressure_val = Tube.State1.p + + # Get the shell pressure based on either the inlet or outlet pressure + # based on whether it is a low-pressure or high-pressure shell + if shell_pressure == 'low': + back_pressure = min((inlet_pressure, outlet_pressure)) + elif shell_pressure == 'low:shell': + back_pressure = min((shell_pressure_val, outlet_pressure)) + elif shell_pressure == 'high': + back_pressure = max((inlet_pressure, outlet_pressure)) + elif shell_pressure == 'mid': + back_pressure = (inlet_pressure + outlet_pressure)/2 + elif shell_pressure == 'mean:shell': + #Average Pressure shell top of vane + back_pressure = np.mean(self.CVs['shell'].State.p*1000) + else: + raise KeyError("keyword argument shell_pressure must be one of 'low', 'low:shell', 'mid' or 'high'; received '"+str(shell_pressure)+"'") + + self.calculate_force_terms(p_shell = back_pressure) + self.Wdot_friction = self.Ler + self.Lrc + self.Lec + self.Lv + self.Ls + self.Lbs + self.Lbl + + Nstep = self.Itheta+1 + _slice = range(self.Itheta+1) + + self.Ler_mean = (np.trapz(self.Ler, self.t[_slice])/(2*pi)) + self.Lrc_mean = (np.trapz(self.Lrc, self.t[_slice])/(2*pi)) + self.Lec_mean =(np.trapz(self.Lec, self.t[_slice])/(2*pi)) + self.Lv_mean = (np.trapz(self.Lv, self.t[_slice])/(2*pi)) + self.Ls_mean = (np.trapz(self.Ls, self.t[_slice])/(2*pi)) + self.Lbs_mean = (np.trapz(self.Lbs, self.t[_slice])/(2*pi)) + self.Lbl_mean = (np.trapz(self.Lbl, self.t[_slice])/(2*pi)) + + self.Wdot_friction_mean = (np.trapz(self.Wdot_friction, self.t[_slice])/(2*pi))/1000 + print('average mechanical and friction losses [kW]:',self.Wdot_friction_mean) + + return self.Wdot_friction_mean + + def mechanical_losses(self): + + if hasattr(self.mech,'detailed_analysis') and self.mech.detailed_analysis == True: + #Friction losses + Wdot_friction = self.friction_losses() #[kW] + return Wdot_friction + else: + print('Simplified mechanical losses') + return self.mech.Wdot_parasitic #[kW] + + def ambient_heat_transfer(self): + """ + The ambient heat transfer for the compressor in kW + + Returns a positive value if heat is added to the compressor from the + ambient + + The heat transfer is obtained for natural convection from a vertical cylinder + """ + ref = 'Air' + T_w = self.Tlumps[0] # [K] + T = self.Tamb # [K] + P_film = 101.325 #[kPa] + + # Vertical side shell + h_c_vertical = HTC('vertical_plate',T_w,T,P_film,ref,self.Hshell) + A_shell_vertical = 2*pi*self.Rshell*self.Hshell # m2 + Q_vertical = h_c_vertical*A_shell_vertical*(T-T_w) + + # Top Surface + h_c_topshell = HTC('horizontal_plate',T_w,T,P_film,ref,self.Rshell*2,PlateNum='upper_surface') + A_shell_top = pi*self.Rshell**2 # m2 + Q_topshell = h_c_topshell*A_shell_top*(T-T_w) + + # Bottom Surface + h_c_bottomshell = HTC('horizontal_plate',T_w,T,P_film,ref,self.Rshell*2,PlateNum='lower_surface') + A_shell_bottom = pi*self.Rshell**2 # m2 + Q_bottomshell = h_c_bottomshell*A_shell_bottom*(T-T_w) + + # Total heat transfer + Q = Q_vertical + Q_topshell + Q_bottomshell # [kW] + + return Q + + def lump_energy_balance_callback(self): + """ + A callback used in PDSimCore.solve to do the energy balance on the lump + + Note: we neglect heat transfer to the gas in the working chamber + """ + #Mechanical and friction losses are added to the lump + self.Wdot_losses = self.mechanical_losses() #[kW] + + #Heat transfer tubes with shell + self.Qtubes = -sum([Tube.Q for Tube in self.Tubes]) + + #Heat transfer between the shell and the ambient + self.Qamb = self.ambient_heat_transfer() #[kW] + + #Mechanical power + self.Wdot_mechanical = self.Wdot_pv + self.Wdot_losses + + #The actual torque required to do the compression [N-m] + self.tau_mechanical = self.Wdot_mechanical / self.omega * 1000 + + #Mechanical efficiency + self.eta_mechanical = self.Wdot_pv/self.Wdot_mechanical + + if self.motor.type == 'const_eta_motor': + self.eta_motor = self.motor.eta_motor + elif self.motor.type == 'motor_map': + # Use the motor map to calculate the slip rotational speed [rad/s] + # and the motor efficiency as a function of the torque [N-m] + eta, omega = self.motor.apply_map(self.tau_mechanical) + self.eta_motor = eta + self.omega = omega + else: + raise AttributeError('motor.type must be one of "const_eta_motor" or "motor_map"') + + #Motor losses [kW] + self.motor.losses = self.Wdot_mechanical*(1/self.eta_motor-1) + + #Electrical Power + self.Wdot_electrical = self.Wdot_mechanical + self.motor.losses + + if hasattr(self,'Wdot_i'): + #Overall isentropic efficiency + self.eta_oi = self.Wdot_i/self.Wdot_electrical + + if self.verbosity > 0: + print('At this iteration') + print(' Indicated power:', self.Wdot_pv,'kW') + print(' Mechanical torque:',self.tau_mechanical,'Nm') + print(' Electrical power:', self.Wdot_electrical,'kW') + print(' Mechanical losses:', self.Wdot_losses,'kW') + print(' Motor losses:', self.motor.losses,'kW') + print(' Mass flow rate:', self.mdot,'kg/s', self.mdot*3600,'kg/h') + if hasattr(self,'Wdot_i'): + print(' Over. isentropic:', self.eta_oi,'-') + if hasattr(self,'eta_v'): + print(' Volumetric:', self.eta_v,'-') + if self.mech.detailed_analysis: + print('Mechanical efficiency:',self.eta_mechanical,'-') + print('Ecc. & inner roller :', self.Ler_mean,'W') + print('Roller face & cylinder head :', self.Lrc_mean,'W') + print('Ecc. face & cylinder head :', self.Lec_mean,'W') + print('Vane tip & roller :', self.Lv_mean,'W') + print('Vane sides & vane slot :', self.Ls_mean,'W') + print('Shaft & upper bearing :', self.Lbs_mean,'W') + print('Shaft & lower bearing :', self.Lbl_mean,'W') + + return self.Qamb + self.Qtubes + 0.1*self.motor.losses + + def step_callback(self,t,h,Itheta): + """ + Here we test whether the control volumes need to be + a) Merged + b) Adjusted because you are at the discharge angle + + """ + if hasattr(self.mech,'cylinder_config') and self.mech.cylinder_config == 'dual-cylinder': + if self.solver == 'RK45': + + # This gets called at every step, or partial step + self.theta = t + + def angle_difference(angle1,angle2): + # Due to the periodicity of angles, you need to handle the case where the + # angles wrap around - suppose theta_d is 6.28 and you are at an angles of 0.1 rad + #, the difference should be around 0.1, not -6.27 + # + # This brilliant method is from http://blog.lexique-du-net.com/index.php?post/Calculate-the-real- + #difference-between-two-angles-keeping-the-sign + # and the comment of user tk + return (angle1-angle2+pi)%(2*pi)-pi + + disable=False + + if t pi: + #print 'here h:',h,t + disable = False + # if t > 3.5: + # raise ValueError + return disable,h + else: + self.theta = t + return False, h + else: + self.theta = t + return False, h + + + + diff --git a/examples/rollingpiston_compressor.py b/examples/rollingpiston_compressor.py new file mode 100644 index 0000000..2290036 --- /dev/null +++ b/examples/rollingpiston_compressor.py @@ -0,0 +1,263 @@ +""" +Domestic Hermetic Rolling Piston Compressor + +References: +K.T. Ooi, "Design optimization of a rolling piston compressor for refrigerators", Applied Thermal Engineering, 25(2005), 813-829 + +K.T. Ooi, "Compressor Performance Comparison When Using R134a nad R1234yf as Working Fluids", International Compressor Engineering Conference at Purdue, 2012. paper 2146 +""" +## This file is meant for experimentation with PDSim features and may not +## work for you. Caveat emptor!! +# +from __future__ import division, print_function + +import sys, os +from math import pi, cos, sin, asin,tan,sqrt + +#PDSim imports +from PDSim.flow.flow import FlowPath +from PDSim.flow import flow_models +from PDSim.misc.datatypes import arraym,listm +from PDSim.core.containers import ControlVolume, Tube +from PDSim.core.core import PDSimCore,struct +from PDSim.plot.plots import debug_plots +#from PDSim.flow.flow_models import ValveModel +from PDSim.flow.ValveModels import ValveModel +from PDSim.core.motor import Motor +from PDSim.rolling.core_clean import RollingPiston + +#CoolProp imports +from CoolProp import State,AbstractState +from CoolProp import CoolProp as CP + +import numpy as np +import matplotlib.pyplot as plt +import timeit + +def RollingPistonCompressor(Te = -10.6, Tc = 54.4, Ref = 'HEOS::R1234yf',HDF5file='rollingpiston_compressor.h5'): + + rolling=RollingPiston() #Instantiate the model + #This runs if the module code is run directly + + rolling.e = 4.8e-3 #Eccentricity, m + rolling.Rc = 29.0e-3 #Radius of cylinder, m + rolling.Rc_outer = 35.0e-3 #Radius of outer cylinder, m + rolling.Rr = 23.4e-3 #Radius roller, m + rolling.Rr_i = 20e-3 #Inner radius roller, m + rolling.Re = 18.0e-3 #Eccentric inner radius, m + rolling.Rs = 10e-3 #Shaft radius, m + rolling.rv = 2.6e-3 #Radius of vane tip, m + rolling.Hc = 44.0e-3 #Compressor cylinder height, m + rolling.b = 4.7e-3 #Vane thickness, m + rolling.hv = 27.8e-3 #Vane height, m + rolling.L_slot = 15.2e-3 #Length of vane slot, m + rolling.Hshell = 100e-3 #Height of compressor shell, m + rolling.Rshell = 40e-3 #Compressor shell outer diameter, m + rolling.Nmot = 2875 #Electric motor rotational speed, rpm + rolling.omega = 2*pi*rolling.Nmot/60 #Frequency, rad/sec (60Hz) + rolling.d_discharge= 8e-3 #Discharge port diameter, m + rolling.d_suction= 14e-3 #Suction diameter, m + rolling.A_discharge=pi*rolling.d_discharge**2/4 + rolling.A_suction=pi*rolling.d_suction**2/4 + rolling.D_sp = 12e-3 #Suction pipe diameter, m + rolling.D_dpc = 7.5e-3 #Diameter of processing tool for discharge cut, m + + # Half angle of suction port, rad + rolling.delta_theta_s = rolling.D_sp/(2*rolling.Rc) + # Center angle of suction port in cylinder wrt vane, rad + rolling.theta_sm = 22*pi/180 + # Starting angle of suction port in cylinder wrt vane, rad + rolling.theta_s1 = rolling.theta_sm - rolling.delta_theta_s + # Ending angle of suction port in cylinder wrt vane, rad + rolling.theta_s2 = rolling.theta_sm + rolling.delta_theta_s + + # Half angle of discharge port, rad + rolling.delta_theta_d = rolling.D_dpc/(2*rolling.Rc) + # Center angle of discharge port in cylinder wrt vane, rad + rolling.theta_dm_ang = 20*pi/180 + # Center angle of discharge port in cylinder wrt vane, rad + rolling.theta_dm = 2*pi - rolling.theta_dm_ang + # Starting angle of discharge port in cylinder wrt vane, rad + rolling.theta_d1 = rolling.theta_dm - rolling.delta_theta_d + # Ending angle of discharge port in cylinder wrt vane, rad + rolling.theta_d2 = rolling.theta_dm + rolling.delta_theta_d + + #Vertical height from base of tool angle to cylinder end, m + rolling.Rp = 29e-3 #Distance between cylinder center and end of tool angle, m + rolling.theta_p = 59*pi/180 #Angle of processing tool discharge cut, radians + rolling.hp = (rolling.Rp - rolling.Rc)*tan(rolling.theta_p) + + # Leakage gaps + rolling.delta_gap = 50e-6 #[m] + rolling.delta_min = 55e-6 #Counter set clearance (at 270 degrees), m + rolling.delta_max = 25e-6 #Set clearance (at 90 degrees), m + rolling.delta_side = 20e-6 # Vertical clearance between roller and cylinder (one side), m + rolling.delta_vb = 20e-6 # Vertical clearance between vane and cylinder (one side) [m] + rolling.delta_vt = 3e-6 #Radial clearance between vane tip and roller, m + rolling.A_ls = 692.78e-6 #Area from lower shell to upper shell, m^2 + rolling.xs = 0 #[-] + rolling.xc = 0 #[-] + rolling.x_slot = 0.5 #[-] + rolling.x_roller = 0.8 #High concentration of oil + + # Flow coeff + rolling.Cflow_roller = 0.5 #[-] + rolling.Cflow_vb = 0.5 #[-] + rolling.Cflow_vt = 0.5 #[-] + rolling.Cflow_32 = 0.43 #[-] + + #Selection of HT model + rolling.HT_on = True + #These are parameters needed for the ambient heat transfer model + rolling.A_shell = pi*rolling.Hshell*2*rolling.Rshell #[m2] + rolling.Tamb = 25 + 273.15 #[K] + + #Parameters for the mechanical losses model (simplified) + rolling.mech = struct() + rolling.mech.detailed_analysis = True + rolling.mech.Wdot_parasitic = 0.01 #Parasitic losses [kW] + + #Motor + rolling.motor = Motor() + rolling.motor.set_eta(0.8) # Ooi, 2005 + rolling.motor.suction_fraction = 1.0 #[-] + + #Boundary condition + rolling.Ref=Ref + rolling.rho_oil = 820 #[kg/m3] + Te_K = Te + 273.15 #[K] + Tc_K = Tc + 273.15 + DT_sh = 11.1 #[K] + Tin = Te_K + DT_sh #[K] + DT_sc = 7 #[K] + temp = State.State(Ref,{'T':Te_K,'Q':1}) #[K] + pe = temp.p #[kPa] + temp.update(dict(T=Tc_K, Q=1)) + pc = temp.p #[kPa] + inletState = State.State(Ref,{'T':Tin,'P':pe}) + + T2s = rolling.guess_outlet_temp(inletState,pc) #[K] + outletState = State.State(Ref,{'T':T2s,'P':pc}) + + #Guess mass flow rate + rolling.Vdisp = rolling.V_disp(rolling.theta_sm,rolling.theta_dm_ang) #Displacement volume + print('Compressor displacement (cm3/rev):',rolling.Vdisp*1e6) + + mdot_guess = inletState.rho*rolling.Vdisp*rolling.Nmot/60 #[kg/s] + + #Add the inlet tube + rolling.add_tube( Tube(key1='inlet.1', + key2='inlet.2', + L=0.03,ID=14e-3, + mdot=mdot_guess, + State1=inletState.copy(), + fixed=1, + TubeFcn=rolling.TubeCode) + ) + + #Add the outlet tube + rolling.add_tube( Tube(key1='outlet.1', + key2='outlet.2', + L=0.03,ID=8e-3, + mdot=mdot_guess, + State2=outletState.copy(), + fixed=2, + TubeFcn=rolling.TubeCode) + ) + + #Add the control volumes. + """ + 'A' = Suction CV + 'B' = Compression/discharge CV + 'shell' = compressor shell + """ + rolling.add_CV( ControlVolume(key='A', + initialState=inletState.copy(), + VdVFcn=rolling.Vs_dVs_ParkYC, + becomes='B') ) + rolling.add_CV( ControlVolume(key='B', + initialState=inletState.copy(), + VdVFcn=rolling.Vc_dVc_ParkYC, + becomes='A') ) + + rolling.add_CV( ControlVolume(key='shell', + initialState=outletState.copy(), + VdVFcn=rolling.V_shell, + becomes='shell') ) + + #Add the flow paths that link flow nodes together + rolling.add_flow(FlowPath(key1='inlet.2',key2='A',MdotFcn=rolling.Suction)) + rolling.add_flow(FlowPath(key1='B',key2='shell',MdotFcn=rolling.Discharge)) + rolling.add_flow(FlowPath(key1='shell',key2='outlet.1',MdotFcn=rolling.Outlet)) + rolling.add_flow(FlowPath(key1='A',key2='B',MdotFcn=rolling.CompressibleGasLeakage_vb)) + rolling.add_flow(FlowPath(key1='A',key2='B',MdotFcn=rolling.CompressibleGasLeakage_vt)) + rolling.add_flow(FlowPath(key1='A',key2='B',MdotFcn=rolling.CompressibleGasLeakage_32)) + + #Add the discharge valve parameters + E = 1.93e11 #Youngs Modulus, [Pa] + h_valve = 0.2e-3 #Valve thickness, [m] + rho_valve = 7200 #Density of spring steel, [kg/m^3] + C_D = 1.17 #Drag coefficient, [-] + d_valve_discharge = 10e-3 #Discharge valve diameter, [m] + l_valve_discharge = 21e-3 #Total valve length, [m] + a_valve_discharge = 16.5e-3 #Distance from anchor point to center of discharge port, [m] + x_stopper_discharge = 3.0e-3 #Stopper height at center of discharge port, [m] + + I=(d_valve_discharge*h_valve**3)/12 #Moment of Intertia for discharge valve,[m^4] + k_valve= (6*E*I)/(a_valve_discharge**2*(3*l_valve_discharge-a_valve_discharge)) #Valve stiffness + m_eff=(1/3)*rho_valve*l_valve_discharge*d_valve_discharge*h_valve #Effective mass of valve reeds, [kg] + x_tr_discharge = 0.25*(rolling.d_discharge**2/d_valve_discharge) #Critical lift position, [m] + + #Define discharge valve + rolling.discharge_valve=ValveModel( + d_valve=d_valve_discharge, + d_port=rolling.d_discharge, + C_D= 1.14, + rho_valve=rho_valve, + x_stopper=x_stopper_discharge, + m_eff = m_eff, + k_valve = k_valve, + x_tr = x_tr_discharge, + key_up='B', + key_down='shell' + ) + #Add the discharge valve + rolling.add_valve(rolling.discharge_valve) + + #Connect the callbacks for the step, endcycle, heat transfer and lump energy balance + rolling.connect_callbacks(step_callback = rolling.step_callback, + endcycle_callback=rolling.endcycle_callback, + heat_transfer_callback=rolling.heat_transfer_callback, + lumps_energy_balance_callback = rolling.lump_energy_balance_callback + ) + #Set debug verbosity level + rolling.verbosity = 1 + + t1 = timeit.default_timer() + rolling.solver = 'RK45' + rolling.solve(key_inlet='inlet.1', + key_outlet='outlet.2', + solver_method = rolling.solver, + eps_cycle = 0.01,#0.001, + eps_energy_balance = 0.01, + hmin=1e-5, + max_number_of_steps = 100000, + plot_every_cycle = False + ) + + print('time taken', timeit.default_timer()-t1) + + debug_plots(rolling) + + del rolling.FlowStorage + from PDSim.misc.hdf5 import HDF5Writer + h5 = HDF5Writer() + h5.write_to_file(rolling, HDF5file) + + return rolling + + +if __name__=='__main__': + + RollingPistonCompressor() + \ No newline at end of file