Merge pull request #78 from ibell/HEOSmix

HEOS mix

PDSim/core/core.py

@@ -24,6 +24,8 @@
     from scipy.integrate import trapz
 except ImportError:
     from PDSim.misc.scipylike import trapz
+
+import scipy.optimize
 
 import h5py
 
@@ -456,10 +458,44 @@ def guess_outlet_temp(self, inlet_state, p_outlet, eta_a=0.7):
         This function can also be overloaded by the subclass in order to 
         implement a different guess method
         """
-
-        h1 = inlet_state.h
-        out_state = inlet_state.copy()
-        out_state.update(dict(S = inlet_state.s, P = p_outlet))
+        try:
+            # For an adiabatic process in which the ratio of heat capacities 
+            # gamma is approximately constant, you can calculate the 
+            # relation between T and p along the isentropic path from:
+            # p1^(1-gamma)*T1^gamma = p2^(1-gamma)*T2^gamma
+
+            h1 = inlet_state.h
+            s1 = inlet_state.s 
+            out_state = inlet_state.copy()
+
+            # See: https://en.wikipedia.org/wiki/Table_of_thermodynamic_equations#Thermodynamic_processes
+            gamma = inlet_state.cp/inlet_state.cv
+            T2s = (inlet_state.p**(1-gamma)*inlet_state.T**gamma/p_outlet**(1-gamma))**(1/gamma)
+            # Solver to correct for non-constant gamma to get isentropic enthalpy h2s, 
+            # should be a small correction
+            def objective_h2s(T):
+                out_state.update(dict(T=T, P=p_outlet))
+                return out_state.s - s1
+            T2s = scipy.optimize.newton(objective_h2s, T2s)
+            h2s = out_state.h
+
+            if p_outlet > inlet_state.p:
+                # Compressor Mode
+                h2 = h1 + (h2s-h1)/eta_a
+            else:
+                # Expander Mode
+                h2 = h1 + (h2s-h1)*eta_a
+            # Iterate with Newton's method to get temperature corresponding
+            # to specified h2
+            def objective_h2(T):
+                out_state.update(dict(T=T, P=p_outlet))
+                return out_state.h - h2
+            return scipy.optimize.newton(objective_h2, T2s)
+
+        except BaseException as be:
+            h1 = inlet_state.h
+            out_state = inlet_state.copy()
+            out_state.update(dict(S = inlet_state.s, P = p_outlet))
         h2s = out_state.h
         if p_outlet > inlet_state.p:
             # Compressor Mode
@@ -748,7 +784,11 @@ def post_cycle(self):
         # for some conditions and you are better off just calculating the enthalpy
         # directly
         temp = outletState.copy()
-        temp.update(dict(P=outletState.p, S=s1))
+        def objective(T):
+            temp.update(dict(P=outletState.p, T=T))
+            return temp.s - s1
+        scipy.optimize.newton(objective, outletState.T)
+
         h2s = temp.h
 
         if outletState.p > inletState.p:
@@ -1211,7 +1251,7 @@ def OBJECTIVE_CYCLE(self, Td_Tlumps0, X, epsilon_cycle = 0.003, epsilon_energy_b
 
                 if len(self.solvers.hdisc_history) == 1:
                     # The first time we get here, perturb the discharge enthalpy
-                    self.Tubes.Nodes[self.key_outlet].update_ph(self.Tubes.Nodes[self.key_outlet].p, h_outlet + 5)
+                    h_target = h_outlet + 5
                 else:
                     #  Get the values from the history
                     hdn1, EBn1 = self.solvers.hdisc_history[-1]
@@ -1222,7 +1262,13 @@ def OBJECTIVE_CYCLE(self, Td_Tlumps0, X, epsilon_cycle = 0.003, epsilon_energy_b
 
                     #  Reset the outlet enthalpy of the outlet tube based on our new
                     #  value for it
-                    self.Tubes.Nodes[self.key_outlet].update_ph(self.Tubes.Nodes[self.key_outlet].p, hdnew)
+                    h_target = hdnew
+
+                # Iterate to solve the H, P flash calculation
+                def objective(T):
+                    self.Tubes.Nodes[self.key_outlet].update(dict(T=T, P=self.Tubes.Nodes[self.key_outlet].p))
+                    return self.Tubes.Nodes[self.key_outlet].h - h_target
+                scipy.optimize.newton(objective, self.Tubes.Nodes[self.key_outlet].T)
 
                 print(self.solvers.hdisc_history)
 

PDSim/flow/flow_models.pyx

@@ -3,6 +3,8 @@ from __future__ import division
 import cython
 cimport cython
 
+from CoolProp cimport constants_header
+
 #Uncomment this line to use the python math functions
 #from math import log,pi,e,pow,sqrt
 
@@ -264,11 +266,35 @@ cpdef IsothermalWallTube(mdot,State1,State2,fixed,L,ID,OD=None,HTModel='Twall',T
             T2_star = T_wall-(T_wall-T1)*exp(-pi*ID*L*alpha/(mdot*cp))
 
             # Get the actual outlet enthalpy based on the additional heat input
-            S_star = State(Fluid,{'T':T2_star,'P':p + DELTAP/1000.0})
+            S_star = State1.copy()
+            S_star.update({'T':T2_star,'P':p + DELTAP/1000.0})
 
             h2 = S_star.h + Q_add/mdot/1000.0
 
-            State2.update({'H':h2,'P':p+DELTAP/1000})
+            try:
+                # The normal P, H calculation, will work for REFPROP and pure and
+                # pseudo-pure fluids in CoolProp
+                State2.update({'H':h2, 'P':p+DELTAP/1000})
+            except:
+                # Special-case mixtures for the calculation of H,P inputs
+                # because CoolProp basically doesn't support H,P inputs
+                # so a local Newton's method approach is used instead
+                # and we cannot use scipy.optimize.newton because Cython
+                # does not support closures in cpdef functions.
+                if hasattr(State2,'pAS') and len(State2.pAS.fluid_names()) > 1:
+                    # Use Newton's method to solve for the 
+                    # temperature yielding the outlet enthalpy
+                    T = T2_star
+                    for counter in range(30):
+                        State2.update(dict(T=T, P=p+DELTAP/1000))
+                        r = State2.h - h2
+                        rprime = State2.pAS.first_partial_deriv(constants_header.iHmass, constants_header.iT, constants_header.iP)/1000 # AbstractState has enthalpy in kJ/kg
+                        dT = -r/rprime
+                        T += dT
+                        if abs(r) < 0.02:
+                            break
+                else :
+                    raise ValueError("unclear what to do here")
 
             # Q is defined to be positive if heat transferred from wall to fluid
             #

PDSim/scroll/core.py

@@ -18,6 +18,7 @@
 from math import pi, exp
 import numpy as np
 import copy
+import scipy.optimize
 
 # If scipy is available, use its interpolation and optimization functions, otherwise, 
 # use our implementation (for packaging purposes mostly)
@@ -949,7 +950,8 @@ def auto_add_CVs(self,inletState,outletState):
                 #It is the outermost pair of compression chambers
                 initState = State.State(inletState.Fluid,
                                         dict(T=inletState.T,
-                                             D=inletState.rho)
+                                             D=inletState.rho), 
+                                        phase='gas'
                                         )
 
             else:
@@ -1866,17 +1868,35 @@ def initial_motor_losses(self, eta_a = 0.8):
 
         inletState = self.Tubes.Nodes[self.key_inlet]
         outletState = self.Tubes.Nodes[self.key_outlet]
+
+        h1 = inletState.h
         s1 = inletState.s
+        outletState = inletState.copy()
+
+        # For an adiabatic process in which the ratio of heat capacities 
+        # gamma is approximately constant, you can calculate the 
+        # relation between T and p along the isentropic path from:
+        # p1^(1-gamma)*T1^gamma = p2^(1-gamma)*T2^gamma
+
+        # See: https://en.wikipedia.org/wiki/Table_of_thermodynamic_equations#Thermodynamic_processes
+        gamma = inletState.cp/inletState.cv
+        T2s = (inletState.p**(1-gamma)*inletState.T**gamma/outletState.p**(1-gamma))**(1/gamma)
+
+        # Solver to correct for non-constant gamma to get isentropic enthalpy h2s, 
+        # should be a small correction
+        def objective_h2s(T):
+            outletState.update(dict(T=T, P=outletState.p))
+            return outletState.s - s1
+        T2s = scipy.optimize.newton(objective_h2s, T2s)
+        h2s = outletState.h
+
+        # And now we iterate to get the given efficiency
         h1 = inletState.h
-        temp = inletState.copy()
-        temp.update(dict(S = s1, P = outletState.p))
-        h2s = temp.h
-
         if outletState.p > inletState.p:
-            #Compressor Mode
+            # Compressor Mode
             h2 = h1 + (h2s-h1)/eta_a
         else:
-            #Expander Mode
+            # Expander Mode
             h2 = h1 + (h2s-h1)*eta_a
 
         # A guess for the compressor mechanical power based on 70% efficiency [kW]

conda_environment.yml

@@ -12,5 +12,5 @@ dependencies:
   - pip
   - pip:
     - CoolProp
-    - Cython
+    - Cython<3
     - quantities

examples/scroll_compressor.py

@@ -68,14 +68,14 @@ def Compressor(ScrollClass, Te = 273, Tc = 300, f = None, OneCycle = False, Ref
     Te = -20 + 273.15
     Tc = 20 + 273.15
     Tin = Te + 11.1
-    temp = State.State(Ref,{'T':Te,'Q':1})
+    temp = State.State(Ref,{'T':Te,'Q':1}, phase='gas')
     pe = temp.p
     temp.update(dict(T=Tc, Q=1))
     pc = temp.p
-    inletState = State.State(Ref,{'T':Tin,'P':pe})
+    inletState = State.State(Ref,{'T':Tin,'P':pe}, phase='gas')
 
-    T2s = ScrollComp.guess_outlet_temp(inletState,pc)
-    outletState = State.State(Ref,{'T':T2s,'P':pc})
+    T2s = ScrollComp.guess_outlet_temp(inletState, pc)
+    outletState = State.State(Ref,{'T':T2s,'P':pc}, phase='gas')
 
     mdot_guess = inletState.rho*ScrollComp.Vdisp*ScrollComp.omega/(2*pi)
 

examples/scroll_compressor_bicubic_with_mixtures.py

@@ -1,15 +1,21 @@
-from __future__ import print_function
+import os 
 
 from PDSim.scroll.core import Scroll
 from scroll_compressor import Compressor
 
-# from CoolProp import CoolProp as CP
-# CP.set_config_string(CP.ALTERNATIVE_REFPROP_PATH, r'/home/ian/REFPROP911')
+from CoolProp import CoolProp as CP
+CP.set_config_string(CP.ALTERNATIVE_REFPROP_PATH, os.getenv('RPPREFIX'))
 
 # CP.set_debug_level(1)
 
-# Use this one to use R410A as a mixture
-Compressor(Scroll, Ref='BICUBIC&REFPROP::R32[0.697614699375863]&R125[0.302385300624138]')
+# Use this one to use R410A as a pseudo-pure fluid
+Compressor(Scroll, Ref='R410A', HDF5file='R410APPF.h5')
 
-# Use this one to use R410A as a pure fluid
-#Compressor(Ref='R410A')
+# Use this one to use R410A as a mixture w/ REFPROP
+Compressor(Scroll, Ref='REFPROP::R32[0.697614699375863]&R125[0.302385300624138]', HDF5file='REFPROPR410Amix.h5')
+
+# Use this one to use R410A as a mixture w/ REFPROP and bicubic interpolation
+Compressor(Scroll, Ref='BICUBIC&REFPROP::R32[0.697614699375863]&R125[0.302385300624138]', HDF5file='REFPROPBICUBICR410Amix.h5')
+
+# Use this one to use R410A as a mixture w/ the HEOS backend of CoolProp
+Compressor(Scroll, Ref='HEOS::R32[0.697614699375863]&R125[0.302385300624138]', HDF5file='HEOSR410Amix.h5')