search

CalebBell / thermo

Thermodynamics and Phase Equilibrium component of Chemical Engineering Design Library (ChEDL)
MIT License
594 stars 114 forks source link

Receive the following error when trying to calculate temperature from pressure and entropy #11

Closed Pheenixm closed 3 years ago

Pheenixm commented 6 years ago 
Traceback (most recent call last):
  File "/Users/RobMontgomery/Documents/Programming/workspace/ChillSeekers/Ejector.py", line 183, in <module>
    flow1, velocity_outlet = ejector.calculate_nozzle_outlet_flow(flow_throat, velocity_throat)
  File "/Users/RobMontgomery/Documents/Programming/workspace/ChillSeekers/Ejector.py", line 106, in calculate_nozzle_outlet_flow
    isoentropic_enthalpy_outlet = Chemical(self.medium, P=assumed_outlet_pressure, T=Chemical(self.medium).calculate_PS(assumed_outlet_pressure, nozzle_flow_entropy)).H
  File "/Users/RobMontgomery/anaconda3/lib/python3.6/site-packages/thermo/chemical.py", line 1107, in calculate_PS
    return newton(to_solve, self.T)
  File "/Users/RobMontgomery/anaconda3/lib/python3.6/site-packages/scipy/optimize/zeros.py", line 184, in newton
    raise RuntimeError(msg)
RuntimeError: Failed to converge after 50 iterations, value is 226.879055881

The code 

    def calculate_nozzle_outlet_flow(self, nozzle_flow, velocity_throat):
        assumed_outlet_pressure = .2e6

        # Calculate Primary Flow Characteristics
        nozzle_flow_entropy = nozzle_flow.S
        nozzle_flow_enthalpy = nozzle_flow.H

        #Calculate flow rate
        mass_flow_nozzle = nozzle_flow.rho * self.nozzle_throat_area * velocity_throat

        iterator = 0
        primary_outlet_flow = None
        while True:
            # Calculate isoentropic enthalpies
            isoentropic_enthalpy_outlet = Chemical(self.medium, P=assumed_outlet_pressure, T=Chemical(self.medium).calculate_PS(assumed_outlet_pressure, nozzle_flow_entropy)).H
            # Determine real enthropy at nozzle
            real_enthalpy_outlet = nozzle_flow_enthalpy - (self.isoentropic_nozzle_efficiency * (nozzle_flow_enthalpy - isoentropic_enthalpy_outlet))
            # Calculate estimated speed of sound
            calculated_mach_velocity = (2 * (nozzle_flow_enthalpy - real_enthalpy_outlet + (velocity_throat ** 2) / 2)) ** 0.5

            # Calculate speed at outlet of nozzle based on known mass flow rate
            primary_outlet_temp = Chemical(self.medium).calculate_PH(assumed_outlet_pressure, real_enthalpy_outlet)
            primary_outlet_flow = Chemical(self.medium, T=primary_outlet_temp, P=assumed_outlet_pressure)
            pof = primary_outlet_flow

            assumed_outlet_rho = pof.rho
            thermo_velocity = mass_flow_nozzle / (self.nozzle_outlet_area * assumed_outlet_rho)

            print('Outlet thermodynamic speed is ', thermo_velocity)
            print('Outlet calculated speed is ', calculated_mach_velocity)
            if math.fabs(thermo_velocity - calculated_mach_velocity) <= 0.1:
                return primary_outlet_flow,thermo_velocity
                break
            else:
                assumed_outlet_pressure = assumed_outlet_pressure - 1000
                iterator += 1
                print(assumed_outlet_pressure)

Any clue what's going on here? I'm trying to perform an iterative conversion for use in a thermodynamic cycle. Thanks.

CalebBell commented 6 years ago

Hello Rob, Ejectors? Very neat stuff! The loop in Chemical to seek the PS values you specified failed to converge. You might want to make a PS plot for the chemical and check the value you're seeking is able to be obtained.

The method Calculate_PS is a very basic newton wrapper of the TP calculation routine; if it works great, but I never had high hopes for it. In particular, newton likes to guess negative temperatures and that obviously goes poorly. Quality is also not part of the API for chemicals at the moment, so you might be running into issues there if it's changing phase. The mixture API may be more helpful there.

Not knowing the chemical you're working with or the PS values you were seeking, I can't speculate further.

Sincerely, Caleb

Pheenixm commented 6 years ago

Here's my entire file, it's for a class, so I have no compunctions about making it publicly available. 

from thermo import Chemical
from thermo import utils
import numpy as np
import math

class Ejector:
    # Geometric Parameters
    nozzle_throat_area = 8.04e-6
    nozzle_outlet_area = 3.914e-5
    mixing_area = 7.088e-5
    diffuser_outlet_area = 4.524e-4
    # Values from Boiler
    primary_flow_temperature = 45 + 273
    # Is pressure temperature dependent?
    primary_flow_pressure = 0.3e6
    # Values from Evaporator
    secondary_flow_temperature = 10 + 273
    secondary_flow_pressure = 0.5e3
    # Estimates
    isoentropic_nozzle_efficiency = 0.95
    isoentropic_suction_efficiency = 0.85
    # Specify fluid
    medium = 'Ammonia'
    ref = "NIST Refprop 9.1"
    entrainment_ratio = 0.5 #secondary mass flow rate / primary mass flow rate

    '''
    Allows you to set all variables manually
    '''
    def __init__(self, 
                primary_flow_temperature=primary_flow_temperature, 
                primary_flow_pressure=primary_flow_pressure, 
                secondary_flow_temperature=secondary_flow_temperature, 
                secondary_flow_pressure=secondary_flow_pressure, 
                nozzle_throat_area=nozzle_throat_area, 
                nozzle_outlet_area=nozzle_outlet_area, 
                mixing_area=mixing_area, 
                diffuser_outlet_area=diffuser_outlet_area, 
                medium=medium,
                entrainment_ratio=entrainment_ratio
                ):

        self.primary_flow_temperature = primary_flow_temperature
        self.primary_flow_pressure = primary_flow_pressure
        self.secondary_flow_temperature = secondary_flow_temperature
        self.secondary_flow_pressure = secondary_flow_pressure
        self.nozzle_throat_area = nozzle_throat_area
        self.nozzle_outlet_area = nozzle_outlet_area
        self.mixing_area = mixing_area
        self.diffuser_outlet_area = diffuser_outlet_area
        self.medium = medium
        self.entrainment_ratio = entrainment_ratio
    '''
    Takes initial conditions and calculates the nozzle pressure from this
    Returns nozzle flow as a Chemical object 
    '''
    def calculate_nozzle_flow(self): 
        # CALCULATE Pressure in the Nozzle (Pt)
        assumed_nozzle_pressure = 0.3e6
        # Calculate Primary Flow Characteristics
        primary_flow = Chemical(self.medium, T=self.primary_flow_temperature, P=self.primary_flow_pressure)
        primary_flow_entropy = primary_flow.S
        primary_flow_enthalpy = primary_flow.H
        iterator = 0
        nozzle_flow = None
        while True:
            # Calculate isoentropic enthalpies
            isoentropic_enthalpy_throat = Chemical(self.medium, P=assumed_nozzle_pressure, T=Chemical(self.medium).calculate_PS(assumed_nozzle_pressure, primary_flow_entropy)).H
            # Determine real enthropy at nozzle
            real_enthalpy_nozzle = primary_flow_enthalpy - (self.isoentropic_nozzle_efficiency * (primary_flow_enthalpy - isoentropic_enthalpy_throat))
            # Calculate estimated speed of sound
            calculated_mach_velocity = math.fabs((primary_flow_enthalpy - real_enthalpy_nozzle) * 2) ** 0.5
            # Calculate speed of sound at the nozzle, using assumed pressure; will have to iterate over this
            primary_nozzle_temp = Chemical(self.medium).calculate_PH(assumed_nozzle_pressure, real_enthalpy_nozzle)
            nozzle_flow = Chemical(self.medium, T=primary_nozzle_temp, P=assumed_nozzle_pressure)
            pnf = nozzle_flow
            dP_dV = pnf.eos.derivatives_and_departures(pnf.T, pnf.P, pnf.Vm, pnf.eos.b, pnf.eos.delta, pnf.eos.epsilon, pnf.eos.a_alpha, pnf.eos.da_alpha_dT, pnf.eos.d2a_alpha_dT2)[0][1]
            mach_velocity = utils.speed_of_sound(pnf.Vm, dP_dV, pnf.Cp, pnf.Cvg, pnf.MW)

            if math.fabs(mach_velocity - calculated_mach_velocity) <= 0.5:
                print("Nozzle temp is %d, Nozzle pressure is %s, Nozzle entropy is %a" % (primary_nozzle_temp, assumed_nozzle_pressure, pnf.S))
                return (nozzle_flow, mach_velocity)
                break
            else:
                assumed_nozzle_pressure = assumed_nozzle_pressure - 1000
                iterator += 1
                # del pnf, primary_nozzle_fluid

    def calculate_nozzle_outlet_flow(self, nozzle_flow, velocity_throat):
        assumed_outlet_pressure = .2e6

        # Calculate Primary Flow Characteristics
        nozzle_flow_entropy = nozzle_flow.S
        nozzle_flow_enthalpy = nozzle_flow.H

        #Calculate flow rate
        mass_flow_nozzle = nozzle_flow.rho * self.nozzle_throat_area * velocity_throat

        iterator = 0
        primary_outlet_flow = None
        while True:
            # Calculate isoentropic enthalpies
            print("Pressure is %d, Entropy is %s" % (assumed_outlet_pressure,nozzle_flow_entropy))
            isoentropic_enthalpy_outlet = Chemical(self.medium, P=assumed_outlet_pressure, T=Chemical(self.medium).calculate_PS(assumed_outlet_pressure, nozzle_flow_entropy)).H
            # Determine real enthropy at nozzle
            real_enthalpy_outlet = nozzle_flow_enthalpy - (self.isoentropic_nozzle_efficiency * (nozzle_flow_enthalpy - isoentropic_enthalpy_outlet))
            # Calculate estimated speed of sound
            calculated_mach_velocity = (2 * (nozzle_flow_enthalpy - real_enthalpy_outlet + (velocity_throat ** 2) / 2)) ** 0.5

            # Calculate speed at outlet of nozzle based on known mass flow rate
            primary_outlet_temp = Chemical(self.medium).calculate_PH(assumed_outlet_pressure, real_enthalpy_outlet)
            primary_outlet_flow = Chemical(self.medium, T=primary_outlet_temp, P=assumed_outlet_pressure)
            pof = primary_outlet_flow

            assumed_outlet_rho = pof.rho
            thermo_velocity = mass_flow_nozzle / (self.nozzle_outlet_area * assumed_outlet_rho)

            print('Outlet thermodynamic speed is ', thermo_velocity)
            print('Outlet calculated speed is ', calculated_mach_velocity)
            if math.fabs(thermo_velocity - calculated_mach_velocity) <= 0.1:
                return primary_outlet_flow,thermo_velocity
                break
            else:
                assumed_outlet_pressure = assumed_outlet_pressure - 1000
                iterator += 1
                print(assumed_outlet_pressure)

        #Flow1 
    def calculate_secondary_pressure(self, flow_outlet, flow_outlet_velocity):
        secondary_flow = Chemical(self.medium, T=self.secondary_flow_temperature, P=self.secondary_flow_pressure)
        secondary_flow_entropy = secondary_flow.S
        secondary_flow_enthalpy = secondary_flow.H
        assumed_suction_pressure = self.secondary_flow_pressure

        while True:
            #First, the S2 Flow#
            isoentropic_enthalpy_suction = Chemical(self.medium, P=assumed_suction_pressure, T=Chemical(self.medium).calculate_PS(assumed_suction_pressure, secondary_flow_entropy)).H
            real_enthalpy_suction = secondary_flow_enthalpy - (self.isoentropic_suction_efficiency * (secondary_flow_enthalpy - isoentropic_enthalpy_suction))

            suction_temp = Chemical(self.medium).calculate_PH(assumed_suction_pressure, real_enthalpy_suction)
            suction_flow = Chemical(self.medium, T=suction_temp, P=assumed_suction_pressure)
            sf = suction_flow
            dP_dV = sf.eos.derivatives_and_departures(sf.T, sf.P, sf.Vm, sf.eos.b, sf.eos.delta, sf.eos.epsilon, sf.eos.a_alpha, sf.eos.da_alpha_dT, sf.eos.d2a_alpha_dT2)[0][1]
            suction_velocity = utils.speed_of_sound(sf.Vm, dP_dV, sf.Cp, sf.Cvg, sf.MW)

            suction_density= sf.rho
            suction_mass_flow = self.entrainment_ratio * (flow_outlet.rho * self.nozzle_outlet_area * flow_outlet_velocity)
            thermo_suction_area = suction_mass_flow / (suction_velocity * suction_density)

            #Now, onto the P2 flow# 
            isoentropic_enthalpy_primary = Chemical(self.medium, P=assumed_suction_pressure, T=Chemical(self.medium).calculate_PS(assumed_suction_pressure, flow_outlet.S)).H
            real_enthalpy_primary = flow_outlet.H - (self.isoentropic_nozzle_efficiency * (flow_outlet.H - isoentropic_enthalpy_primary))

            primary_temp = Chemical(self.medium).calculate_PH(assumed_suction_pressure, real_enthalpy_primary)
            primary_flow_2 = Chemical(self.medium, T=primary_temp, P=assumed_suction_pressure)
            pf = primary_flow_2
            primary_density = pf.rho 
            thermo_primary_area = self.mixing_area - thermo_suction_area
            primary_mass_flow = flow_outlet.rho * self.nozzle_outlet_area * flow_outlet_velocity
            primary_velocity = primary_mass_flow / (thermo_primary_area * primary_density)

            right_hand = primary_mass_flow * (flow_outlet.H + (flow_outlet_velocity ** 2) * 0.5) + suction_mass_flow * secondary_flow_enthalpy
            left_hand = primary_mass_flow * (pf.H + (primary_velocity ** 2) * 0.5) + suction_mass_flow * (sf.H + (suction_velocity ** 2) * 0.5)

            print("Right hand", right_hand)
            print("Left hand" , left_hand)
            if math.fabs(right_hand - left_hand) <= 1000:
                return pf #Primary flow contains all relevant information
                break
            else:
                assumed_suction_pressure = assumed_suction_pressure - 100
                print('Assumed Pressure ', assumed_suction_pressure)            

ejector = Ejector()
flow_throat, velocity_throat = ejector.calculate_nozzle_flow()
print('Pressure in Nozzle Throat is ', flow_throat.P)

flow1, velocity_outlet = ejector.calculate_nozzle_outlet_flow(flow_throat, velocity_throat)
print('Calling first')

print('Pressure in Nozzle Outlet is ', flow1.P)
print('Mass flow rate is ', flow_throat.rho * ejector.nozzle_throat_area * velocity_throat)
mixing_flow = ejector.calculate_secondary_pressure(flow1, velocity_outlet)

from thermo import Chemical
from thermo import utils
import numpy as np
import math

class Ejector:
    # Geometric Parameters
    nozzle_throat_area = 8.04e-6
    nozzle_outlet_area = 3.914e-5
    mixing_area = 7.088e-5
    diffuser_outlet_area = 4.524e-4
    # Values from Boiler
    primary_flow_temperature = 45 + 273
    # Is pressure temperature dependent?
    primary_flow_pressure = 0.3e6
    # Values from Evaporator
    secondary_flow_temperature = 10 + 273
    secondary_flow_pressure = 0.5e3
    # Estimates
    isoentropic_nozzle_efficiency = 0.95
    isoentropic_suction_efficiency = 0.85
    # Specify fluid
    medium = 'Ammonia'
    ref = "NIST Refprop 9.1"
    entrainment_ratio = 0.5 #secondary mass flow rate / primary mass flow rate

    '''
    Allows you to set all variables manually
    '''
    def __init__(self, 
                primary_flow_temperature=primary_flow_temperature, 
                primary_flow_pressure=primary_flow_pressure, 
                secondary_flow_temperature=secondary_flow_temperature, 
                secondary_flow_pressure=secondary_flow_pressure, 
                nozzle_throat_area=nozzle_throat_area, 
                nozzle_outlet_area=nozzle_outlet_area, 
                mixing_area=mixing_area, 
                diffuser_outlet_area=diffuser_outlet_area, 
                medium=medium,
                entrainment_ratio=entrainment_ratio
                ):

        self.primary_flow_temperature = primary_flow_temperature
        self.primary_flow_pressure = primary_flow_pressure
        self.secondary_flow_temperature = secondary_flow_temperature
        self.secondary_flow_pressure = secondary_flow_pressure
        self.nozzle_throat_area = nozzle_throat_area
        self.nozzle_outlet_area = nozzle_outlet_area
        self.mixing_area = mixing_area
        self.diffuser_outlet_area = diffuser_outlet_area
        self.medium = medium
        self.entrainment_ratio = entrainment_ratio
    '''
    Takes initial conditions and calculates the nozzle pressure from this
    Returns nozzle flow as a Chemical object 
    '''
    def calculate_nozzle_flow(self): 
        # CALCULATE Pressure in the Nozzle (Pt)
        assumed_nozzle_pressure = 0.3e6
        # Calculate Primary Flow Characteristics
        primary_flow = Chemical(self.medium, T=self.primary_flow_temperature, P=self.primary_flow_pressure)
        primary_flow_entropy = primary_flow.S
        primary_flow_enthalpy = primary_flow.H
        iterator = 0
        nozzle_flow = None
        while True:
            # Calculate isoentropic enthalpies
            isoentropic_enthalpy_throat = Chemical(self.medium, P=assumed_nozzle_pressure, T=Chemical(self.medium).calculate_PS(assumed_nozzle_pressure, primary_flow_entropy)).H
            # Determine real enthropy at nozzle
            real_enthalpy_nozzle = primary_flow_enthalpy - (self.isoentropic_nozzle_efficiency * (primary_flow_enthalpy - isoentropic_enthalpy_throat))
            # Calculate estimated speed of sound
            calculated_mach_velocity = math.fabs((primary_flow_enthalpy - real_enthalpy_nozzle) * 2) ** 0.5
            # Calculate speed of sound at the nozzle, using assumed pressure; will have to iterate over this
            primary_nozzle_temp = Chemical(self.medium).calculate_PH(assumed_nozzle_pressure, real_enthalpy_nozzle)
            nozzle_flow = Chemical(self.medium, T=primary_nozzle_temp, P=assumed_nozzle_pressure)
            pnf = nozzle_flow
            dP_dV = pnf.eos.derivatives_and_departures(pnf.T, pnf.P, pnf.Vm, pnf.eos.b, pnf.eos.delta, pnf.eos.epsilon, pnf.eos.a_alpha, pnf.eos.da_alpha_dT, pnf.eos.d2a_alpha_dT2)[0][1]
            mach_velocity = utils.speed_of_sound(pnf.Vm, dP_dV, pnf.Cp, pnf.Cvg, pnf.MW)

            if math.fabs(mach_velocity - calculated_mach_velocity) <= 0.5:
                print("Nozzle temp is %d, Nozzle pressure is %s, Nozzle entropy is %a" % (primary_nozzle_temp, assumed_nozzle_pressure, pnf.S))
                return (nozzle_flow, mach_velocity)
                break
            else:
                assumed_nozzle_pressure = assumed_nozzle_pressure - 1000
                iterator += 1
                # del pnf, primary_nozzle_fluid

    def calculate_nozzle_outlet_flow(self, nozzle_flow, velocity_throat):
        assumed_outlet_pressure = .2e6

        # Calculate Primary Flow Characteristics
        nozzle_flow_entropy = nozzle_flow.S
        nozzle_flow_enthalpy = nozzle_flow.H

        #Calculate flow rate
        mass_flow_nozzle = nozzle_flow.rho * self.nozzle_throat_area * velocity_throat

        iterator = 0
        primary_outlet_flow = None
        while True:
            # Calculate isoentropic enthalpies
            print("Pressure is %d, Entropy is %s" % (assumed_outlet_pressure,nozzle_flow_entropy))
            isoentropic_enthalpy_outlet = Chemical(self.medium, P=assumed_outlet_pressure, T=Chemical(self.medium).calculate_PS(assumed_outlet_pressure, nozzle_flow_entropy)).H
            # Determine real enthropy at nozzle
            real_enthalpy_outlet = nozzle_flow_enthalpy - (self.isoentropic_nozzle_efficiency * (nozzle_flow_enthalpy - isoentropic_enthalpy_outlet))
            # Calculate estimated speed of sound
            calculated_mach_velocity = (2 * (nozzle_flow_enthalpy - real_enthalpy_outlet + (velocity_throat ** 2) / 2)) ** 0.5

            # Calculate speed at outlet of nozzle based on known mass flow rate
            primary_outlet_temp = Chemical(self.medium).calculate_PH(assumed_outlet_pressure, real_enthalpy_outlet)
            primary_outlet_flow = Chemical(self.medium, T=primary_outlet_temp, P=assumed_outlet_pressure)
            pof = primary_outlet_flow

            assumed_outlet_rho = pof.rho
            thermo_velocity = mass_flow_nozzle / (self.nozzle_outlet_area * assumed_outlet_rho)

            print('Outlet thermodynamic speed is ', thermo_velocity)
            print('Outlet calculated speed is ', calculated_mach_velocity)
            if math.fabs(thermo_velocity - calculated_mach_velocity) <= 0.1:
                return primary_outlet_flow,thermo_velocity
                break
            else:
                assumed_outlet_pressure = assumed_outlet_pressure - 1000
                iterator += 1
                print(assumed_outlet_pressure)

        #Flow1 
    def calculate_secondary_pressure(self, flow_outlet, flow_outlet_velocity):
        secondary_flow = Chemical(self.medium, T=self.secondary_flow_temperature, P=self.secondary_flow_pressure)
        secondary_flow_entropy = secondary_flow.S
        secondary_flow_enthalpy = secondary_flow.H
        assumed_suction_pressure = self.secondary_flow_pressure

        while True:
            #First, the S2 Flow#
            isoentropic_enthalpy_suction = Chemical(self.medium, P=assumed_suction_pressure, T=Chemical(self.medium).calculate_PS(assumed_suction_pressure, secondary_flow_entropy)).H
            real_enthalpy_suction = secondary_flow_enthalpy - (self.isoentropic_suction_efficiency * (secondary_flow_enthalpy - isoentropic_enthalpy_suction))

            suction_temp = Chemical(self.medium).calculate_PH(assumed_suction_pressure, real_enthalpy_suction)
            suction_flow = Chemical(self.medium, T=suction_temp, P=assumed_suction_pressure)
            sf = suction_flow
            dP_dV = sf.eos.derivatives_and_departures(sf.T, sf.P, sf.Vm, sf.eos.b, sf.eos.delta, sf.eos.epsilon, sf.eos.a_alpha, sf.eos.da_alpha_dT, sf.eos.d2a_alpha_dT2)[0][1]
            suction_velocity = utils.speed_of_sound(sf.Vm, dP_dV, sf.Cp, sf.Cvg, sf.MW)

            suction_density= sf.rho
            suction_mass_flow = self.entrainment_ratio * (flow_outlet.rho * self.nozzle_outlet_area * flow_outlet_velocity)
            thermo_suction_area = suction_mass_flow / (suction_velocity * suction_density)

            #Now, onto the P2 flow# 
            isoentropic_enthalpy_primary = Chemical(self.medium, P=assumed_suction_pressure, T=Chemical(self.medium).calculate_PS(assumed_suction_pressure, flow_outlet.S)).H
            real_enthalpy_primary = flow_outlet.H - (self.isoentropic_nozzle_efficiency * (flow_outlet.H - isoentropic_enthalpy_primary))

            primary_temp = Chemical(self.medium).calculate_PH(assumed_suction_pressure, real_enthalpy_primary)
            primary_flow_2 = Chemical(self.medium, T=primary_temp, P=assumed_suction_pressure)
            pf = primary_flow_2
            primary_density = pf.rho 
            thermo_primary_area = self.mixing_area - thermo_suction_area
            primary_mass_flow = flow_outlet.rho * self.nozzle_outlet_area * flow_outlet_velocity
            primary_velocity = primary_mass_flow / (thermo_primary_area * primary_density)

            right_hand = primary_mass_flow * (flow_outlet.H + (flow_outlet_velocity ** 2) * 0.5) + suction_mass_flow * secondary_flow_enthalpy
            left_hand = primary_mass_flow * (pf.H + (primary_velocity ** 2) * 0.5) + suction_mass_flow * (sf.H + (suction_velocity ** 2) * 0.5)

            print("Right hand", right_hand)
            print("Left hand" , left_hand)
            if math.fabs(right_hand - left_hand) <= 1000:
                return pf #Primary flow contains all relevant information
                break
            else:
                assumed_suction_pressure = assumed_suction_pressure - 100
                print('Assumed Pressure ', assumed_suction_pressure)            

ejector = Ejector()
flow_throat, velocity_throat = ejector.calculate_nozzle_flow()
print('Pressure in Nozzle Throat is ', flow_throat.P)

flow1, velocity_outlet = ejector.calculate_nozzle_outlet_flow(flow_throat, velocity_throat)
print('Calling first')

print('Pressure in Nozzle Outlet is ', flow1.P)
print('Mass flow rate is ', flow_throat.rho * ejector.nozzle_throat_area * velocity_throat)
mixing_flow = ejector.calculate_secondary_pressure(flow1, velocity_outlet)

from thermo import Chemical
from thermo import utils
import numpy as np
import math

class Ejector:
    # Geometric Parameters
    nozzle_throat_area = 8.04e-6
    nozzle_outlet_area = 3.914e-5
    mixing_area = 7.088e-5
    diffuser_outlet_area = 4.524e-4
    # Values from Boiler
    primary_flow_temperature = 45 + 273
    # Is pressure temperature dependent?
    primary_flow_pressure = 0.3e6
    # Values from Evaporator
    secondary_flow_temperature = 10 + 273
    secondary_flow_pressure = 0.5e3
    # Estimates
    isoentropic_nozzle_efficiency = 0.95
    isoentropic_suction_efficiency = 0.85
    # Specify fluid
    medium = 'Ammonia'
    ref = "NIST Refprop 9.1"
    entrainment_ratio = 0.5 #secondary mass flow rate / primary mass flow rate

    '''
    Allows you to set all variables manually
    '''
    def __init__(self, 
                primary_flow_temperature=primary_flow_temperature, 
                primary_flow_pressure=primary_flow_pressure, 
                secondary_flow_temperature=secondary_flow_temperature, 
                secondary_flow_pressure=secondary_flow_pressure, 
                nozzle_throat_area=nozzle_throat_area, 
                nozzle_outlet_area=nozzle_outlet_area, 
                mixing_area=mixing_area, 
                diffuser_outlet_area=diffuser_outlet_area, 
                medium=medium,
                entrainment_ratio=entrainment_ratio
                ):

        self.primary_flow_temperature = primary_flow_temperature
        self.primary_flow_pressure = primary_flow_pressure
        self.secondary_flow_temperature = secondary_flow_temperature
        self.secondary_flow_pressure = secondary_flow_pressure
        self.nozzle_throat_area = nozzle_throat_area
        self.nozzle_outlet_area = nozzle_outlet_area
        self.mixing_area = mixing_area
        self.diffuser_outlet_area = diffuser_outlet_area
        self.medium = medium
        self.entrainment_ratio = entrainment_ratio
    '''
    Takes initial conditions and calculates the nozzle pressure from this
    Returns nozzle flow as a Chemical object 
    '''
    def calculate_nozzle_flow(self): 
        # CALCULATE Pressure in the Nozzle (Pt)
        assumed_nozzle_pressure = 0.3e6
        # Calculate Primary Flow Characteristics
        primary_flow = Chemical(self.medium, T=self.primary_flow_temperature, P=self.primary_flow_pressure)
        primary_flow_entropy = primary_flow.S
        primary_flow_enthalpy = primary_flow.H
        iterator = 0
        nozzle_flow = None
        while True:
            # Calculate isoentropic enthalpies
            isoentropic_enthalpy_throat = Chemical(self.medium, P=assumed_nozzle_pressure, T=Chemical(self.medium).calculate_PS(assumed_nozzle_pressure, primary_flow_entropy)).H
            # Determine real enthropy at nozzle
            real_enthalpy_nozzle = primary_flow_enthalpy - (self.isoentropic_nozzle_efficiency * (primary_flow_enthalpy - isoentropic_enthalpy_throat))
            # Calculate estimated speed of sound
            calculated_mach_velocity = math.fabs((primary_flow_enthalpy - real_enthalpy_nozzle) * 2) ** 0.5
            # Calculate speed of sound at the nozzle, using assumed pressure; will have to iterate over this
            primary_nozzle_temp = Chemical(self.medium).calculate_PH(assumed_nozzle_pressure, real_enthalpy_nozzle)
            nozzle_flow = Chemical(self.medium, T=primary_nozzle_temp, P=assumed_nozzle_pressure)
            pnf = nozzle_flow
            dP_dV = pnf.eos.derivatives_and_departures(pnf.T, pnf.P, pnf.Vm, pnf.eos.b, pnf.eos.delta, pnf.eos.epsilon, pnf.eos.a_alpha, pnf.eos.da_alpha_dT, pnf.eos.d2a_alpha_dT2)[0][1]
            mach_velocity = utils.speed_of_sound(pnf.Vm, dP_dV, pnf.Cp, pnf.Cvg, pnf.MW)

            if math.fabs(mach_velocity - calculated_mach_velocity) <= 0.5:
                print("Nozzle temp is %d, Nozzle pressure is %s, Nozzle entropy is %a" % (primary_nozzle_temp, assumed_nozzle_pressure, pnf.S))
                return (nozzle_flow, mach_velocity)
                break
            else:
                assumed_nozzle_pressure = assumed_nozzle_pressure - 1000
                iterator += 1
                # del pnf, primary_nozzle_fluid

    def calculate_nozzle_outlet_flow(self, nozzle_flow, velocity_throat):
        assumed_outlet_pressure = .2e6

        # Calculate Primary Flow Characteristics
        nozzle_flow_entropy = nozzle_flow.S
        nozzle_flow_enthalpy = nozzle_flow.H

        #Calculate flow rate
        mass_flow_nozzle = nozzle_flow.rho * self.nozzle_throat_area * velocity_throat

        iterator = 0
        primary_outlet_flow = None
        while True:
            # Calculate isoentropic enthalpies
            print("Pressure is %d, Entropy is %s" % (assumed_outlet_pressure,nozzle_flow_entropy))
            isoentropic_enthalpy_outlet = Chemical(self.medium, P=assumed_outlet_pressure, T=Chemical(self.medium).calculate_PS(assumed_outlet_pressure, nozzle_flow_entropy)).H
            # Determine real enthropy at nozzle
            real_enthalpy_outlet = nozzle_flow_enthalpy - (self.isoentropic_nozzle_efficiency * (nozzle_flow_enthalpy - isoentropic_enthalpy_outlet))
            # Calculate estimated speed of sound
            calculated_mach_velocity = (2 * (nozzle_flow_enthalpy - real_enthalpy_outlet + (velocity_throat ** 2) / 2)) ** 0.5

            # Calculate speed at outlet of nozzle based on known mass flow rate
            primary_outlet_temp = Chemical(self.medium).calculate_PH(assumed_outlet_pressure, real_enthalpy_outlet)
            primary_outlet_flow = Chemical(self.medium, T=primary_outlet_temp, P=assumed_outlet_pressure)
            pof = primary_outlet_flow

            assumed_outlet_rho = pof.rho
            thermo_velocity = mass_flow_nozzle / (self.nozzle_outlet_area * assumed_outlet_rho)

            print('Outlet thermodynamic speed is ', thermo_velocity)
            print('Outlet calculated speed is ', calculated_mach_velocity)
            if math.fabs(thermo_velocity - calculated_mach_velocity) <= 0.1:
                return primary_outlet_flow,thermo_velocity
                break
            else:
                assumed_outlet_pressure = assumed_outlet_pressure - 1000
                iterator += 1
                print(assumed_outlet_pressure)

        #Flow1 
    def calculate_secondary_pressure(self, flow_outlet, flow_outlet_velocity):
        secondary_flow = Chemical(self.medium, T=self.secondary_flow_temperature, P=self.secondary_flow_pressure)
        secondary_flow_entropy = secondary_flow.S
        secondary_flow_enthalpy = secondary_flow.H
        assumed_suction_pressure = self.secondary_flow_pressure

        while True:
            #First, the S2 Flow#
            isoentropic_enthalpy_suction = Chemical(self.medium, P=assumed_suction_pressure, T=Chemical(self.medium).calculate_PS(assumed_suction_pressure, secondary_flow_entropy)).H
            real_enthalpy_suction = secondary_flow_enthalpy - (self.isoentropic_suction_efficiency * (secondary_flow_enthalpy - isoentropic_enthalpy_suction))

            suction_temp = Chemical(self.medium).calculate_PH(assumed_suction_pressure, real_enthalpy_suction)
            suction_flow = Chemical(self.medium, T=suction_temp, P=assumed_suction_pressure)
            sf = suction_flow
            dP_dV = sf.eos.derivatives_and_departures(sf.T, sf.P, sf.Vm, sf.eos.b, sf.eos.delta, sf.eos.epsilon, sf.eos.a_alpha, sf.eos.da_alpha_dT, sf.eos.d2a_alpha_dT2)[0][1]
            suction_velocity = utils.speed_of_sound(sf.Vm, dP_dV, sf.Cp, sf.Cvg, sf.MW)

            suction_density= sf.rho
            suction_mass_flow = self.entrainment_ratio * (flow_outlet.rho * self.nozzle_outlet_area * flow_outlet_velocity)
            thermo_suction_area = suction_mass_flow / (suction_velocity * suction_density)

            #Now, onto the P2 flow# 
            isoentropic_enthalpy_primary = Chemical(self.medium, P=assumed_suction_pressure, T=Chemical(self.medium).calculate_PS(assumed_suction_pressure, flow_outlet.S)).H
            real_enthalpy_primary = flow_outlet.H - (self.isoentropic_nozzle_efficiency * (flow_outlet.H - isoentropic_enthalpy_primary))

            primary_temp = Chemical(self.medium).calculate_PH(assumed_suction_pressure, real_enthalpy_primary)
            primary_flow_2 = Chemical(self.medium, T=primary_temp, P=assumed_suction_pressure)
            pf = primary_flow_2
            primary_density = pf.rho 
            thermo_primary_area = self.mixing_area - thermo_suction_area
            primary_mass_flow = flow_outlet.rho * self.nozzle_outlet_area * flow_outlet_velocity
            primary_velocity = primary_mass_flow / (thermo_primary_area * primary_density)

            right_hand = primary_mass_flow * (flow_outlet.H + (flow_outlet_velocity ** 2) * 0.5) + suction_mass_flow * secondary_flow_enthalpy
            left_hand = primary_mass_flow * (pf.H + (primary_velocity ** 2) * 0.5) + suction_mass_flow * (sf.H + (suction_velocity ** 2) * 0.5)

            print("Right hand", right_hand)
            print("Left hand" , left_hand)
            if math.fabs(right_hand - left_hand) <= 1000:
                return pf #Primary flow contains all relevant information
                break
            else:
                assumed_suction_pressure = assumed_suction_pressure - 100
                print('Assumed Pressure ', assumed_suction_pressure)            

ejector = Ejector()
flow_throat, velocity_throat = ejector.calculate_nozzle_flow()
print('Pressure in Nozzle Throat is ', flow_throat.P)

flow1, velocity_outlet = ejector.calculate_nozzle_outlet_flow(flow_throat, velocity_throat)
print('Calling first')

print('Pressure in Nozzle Outlet is ', flow1.P)
print('Mass flow rate is ', flow_throat.rho * ejector.nozzle_throat_area * velocity_throat)
mixing_flow = ejector.calculate_secondary_pressure(flow1, velocity_outlet)

from thermo import Chemical
from thermo import utils
import numpy as np
import math

class Ejector:
    # Geometric Parameters
    nozzle_throat_area = 8.04e-6
    nozzle_outlet_area = 3.914e-5
    mixing_area = 7.088e-5
    diffuser_outlet_area = 4.524e-4
    # Values from Boiler
    primary_flow_temperature = 45 + 273
    # Is pressure temperature dependent?
    primary_flow_pressure = 0.3e6
    # Values from Evaporator
    secondary_flow_temperature = 10 + 273
    secondary_flow_pressure = 0.5e3
    # Estimates
    isoentropic_nozzle_efficiency = 0.95
    isoentropic_suction_efficiency = 0.85
    # Specify fluid
    medium = 'Ammonia'
    ref = "NIST Refprop 9.1"
    entrainment_ratio = 0.5 #secondary mass flow rate / primary mass flow rate

    '''
    Allows you to set all variables manually
    '''
    def __init__(self, 
                primary_flow_temperature=primary_flow_temperature, 
                primary_flow_pressure=primary_flow_pressure, 
                secondary_flow_temperature=secondary_flow_temperature, 
                secondary_flow_pressure=secondary_flow_pressure, 
                nozzle_throat_area=nozzle_throat_area, 
                nozzle_outlet_area=nozzle_outlet_area, 
                mixing_area=mixing_area, 
                diffuser_outlet_area=diffuser_outlet_area, 
                medium=medium,
                entrainment_ratio=entrainment_ratio
                ):

        self.primary_flow_temperature = primary_flow_temperature
        self.primary_flow_pressure = primary_flow_pressure
        self.secondary_flow_temperature = secondary_flow_temperature
        self.secondary_flow_pressure = secondary_flow_pressure
        self.nozzle_throat_area = nozzle_throat_area
        self.nozzle_outlet_area = nozzle_outlet_area
        self.mixing_area = mixing_area
        self.diffuser_outlet_area = diffuser_outlet_area
        self.medium = medium
        self.entrainment_ratio = entrainment_ratio
    '''
    Takes initial conditions and calculates the nozzle pressure from this
    Returns nozzle flow as a Chemical object 
    '''
    def calculate_nozzle_flow(self): 
        # CALCULATE Pressure in the Nozzle (Pt)
        assumed_nozzle_pressure = 0.3e6
        # Calculate Primary Flow Characteristics
        primary_flow = Chemical(self.medium, T=self.primary_flow_temperature, P=self.primary_flow_pressure)
        primary_flow_entropy = primary_flow.S
        primary_flow_enthalpy = primary_flow.H
        iterator = 0
        nozzle_flow = None
        while True:
            # Calculate isoentropic enthalpies
            isoentropic_enthalpy_throat = Chemical(self.medium, P=assumed_nozzle_pressure, T=Chemical(self.medium).calculate_PS(assumed_nozzle_pressure, primary_flow_entropy)).H
            # Determine real enthropy at nozzle
            real_enthalpy_nozzle = primary_flow_enthalpy - (self.isoentropic_nozzle_efficiency * (primary_flow_enthalpy - isoentropic_enthalpy_throat))
            # Calculate estimated speed of sound
            calculated_mach_velocity = math.fabs((primary_flow_enthalpy - real_enthalpy_nozzle) * 2) ** 0.5
            # Calculate speed of sound at the nozzle, using assumed pressure; will have to iterate over this
            primary_nozzle_temp = Chemical(self.medium).calculate_PH(assumed_nozzle_pressure, real_enthalpy_nozzle)
            nozzle_flow = Chemical(self.medium, T=primary_nozzle_temp, P=assumed_nozzle_pressure)
            pnf = nozzle_flow
            dP_dV = pnf.eos.derivatives_and_departures(pnf.T, pnf.P, pnf.Vm, pnf.eos.b, pnf.eos.delta, pnf.eos.epsilon, pnf.eos.a_alpha, pnf.eos.da_alpha_dT, pnf.eos.d2a_alpha_dT2)[0][1]
            mach_velocity = utils.speed_of_sound(pnf.Vm, dP_dV, pnf.Cp, pnf.Cvg, pnf.MW)

            if math.fabs(mach_velocity - calculated_mach_velocity) <= 0.5:
                print("Nozzle temp is %d, Nozzle pressure is %s, Nozzle entropy is %a" % (primary_nozzle_temp, assumed_nozzle_pressure, pnf.S))
                return (nozzle_flow, mach_velocity)
                break
            else:
                assumed_nozzle_pressure = assumed_nozzle_pressure - 1000
                iterator += 1
                # del pnf, primary_nozzle_fluid

    def calculate_nozzle_outlet_flow(self, nozzle_flow, velocity_throat):
        assumed_outlet_pressure = .2e6

        # Calculate Primary Flow Characteristics
        nozzle_flow_entropy = nozzle_flow.S
        nozzle_flow_enthalpy = nozzle_flow.H

        #Calculate flow rate
        mass_flow_nozzle = nozzle_flow.rho * self.nozzle_throat_area * velocity_throat

        iterator = 0
        primary_outlet_flow = None
        while True:
            # Calculate isoentropic enthalpies
            print("Pressure is %d, Entropy is %s" % (assumed_outlet_pressure,nozzle_flow_entropy))
            isoentropic_enthalpy_outlet = Chemical(self.medium, P=assumed_outlet_pressure, T=Chemical(self.medium).calculate_PS(assumed_outlet_pressure, nozzle_flow_entropy)).H
            # Determine real enthropy at nozzle
            real_enthalpy_outlet = nozzle_flow_enthalpy - (self.isoentropic_nozzle_efficiency * (nozzle_flow_enthalpy - isoentropic_enthalpy_outlet))
            # Calculate estimated speed of sound
            calculated_mach_velocity = (2 * (nozzle_flow_enthalpy - real_enthalpy_outlet + (velocity_throat ** 2) / 2)) ** 0.5

            # Calculate speed at outlet of nozzle based on known mass flow rate
            primary_outlet_temp = Chemical(self.medium).calculate_PH(assumed_outlet_pressure, real_enthalpy_outlet)
            primary_outlet_flow = Chemical(self.medium, T=primary_outlet_temp, P=assumed_outlet_pressure)
            pof = primary_outlet_flow

            assumed_outlet_rho = pof.rho
            thermo_velocity = mass_flow_nozzle / (self.nozzle_outlet_area * assumed_outlet_rho)

            print('Outlet thermodynamic speed is ', thermo_velocity)
            print('Outlet calculated speed is ', calculated_mach_velocity)
            if math.fabs(thermo_velocity - calculated_mach_velocity) <= 0.1:
                return primary_outlet_flow,thermo_velocity
                break
            else:
                assumed_outlet_pressure = assumed_outlet_pressure - 1000
                iterator += 1
                print(assumed_outlet_pressure)

        #Flow1 
    def calculate_secondary_pressure(self, flow_outlet, flow_outlet_velocity):
        secondary_flow = Chemical(self.medium, T=self.secondary_flow_temperature, P=self.secondary_flow_pressure)
        secondary_flow_entropy = secondary_flow.S
        secondary_flow_enthalpy = secondary_flow.H
        assumed_suction_pressure = self.secondary_flow_pressure

        while True:
            #First, the S2 Flow#
            isoentropic_enthalpy_suction = Chemical(self.medium, P=assumed_suction_pressure, T=Chemical(self.medium).calculate_PS(assumed_suction_pressure, secondary_flow_entropy)).H
            real_enthalpy_suction = secondary_flow_enthalpy - (self.isoentropic_suction_efficiency * (secondary_flow_enthalpy - isoentropic_enthalpy_suction))

            suction_temp = Chemical(self.medium).calculate_PH(assumed_suction_pressure, real_enthalpy_suction)
            suction_flow = Chemical(self.medium, T=suction_temp, P=assumed_suction_pressure)
            sf = suction_flow
            dP_dV = sf.eos.derivatives_and_departures(sf.T, sf.P, sf.Vm, sf.eos.b, sf.eos.delta, sf.eos.epsilon, sf.eos.a_alpha, sf.eos.da_alpha_dT, sf.eos.d2a_alpha_dT2)[0][1]
            suction_velocity = utils.speed_of_sound(sf.Vm, dP_dV, sf.Cp, sf.Cvg, sf.MW)

            suction_density= sf.rho
            suction_mass_flow = self.entrainment_ratio * (flow_outlet.rho * self.nozzle_outlet_area * flow_outlet_velocity)
            thermo_suction_area = suction_mass_flow / (suction_velocity * suction_density)

            #Now, onto the P2 flow# 
            isoentropic_enthalpy_primary = Chemical(self.medium, P=assumed_suction_pressure, T=Chemical(self.medium).calculate_PS(assumed_suction_pressure, flow_outlet.S)).H
            real_enthalpy_primary = flow_outlet.H - (self.isoentropic_nozzle_efficiency * (flow_outlet.H - isoentropic_enthalpy_primary))

            primary_temp = Chemical(self.medium).calculate_PH(assumed_suction_pressure, real_enthalpy_primary)
            primary_flow_2 = Chemical(self.medium, T=primary_temp, P=assumed_suction_pressure)
            pf = primary_flow_2
            primary_density = pf.rho 
            thermo_primary_area = self.mixing_area - thermo_suction_area
            primary_mass_flow = flow_outlet.rho * self.nozzle_outlet_area * flow_outlet_velocity
            primary_velocity = primary_mass_flow / (thermo_primary_area * primary_density)

            right_hand = primary_mass_flow * (flow_outlet.H + (flow_outlet_velocity ** 2) * 0.5) + suction_mass_flow * secondary_flow_enthalpy
            left_hand = primary_mass_flow * (pf.H + (primary_velocity ** 2) * 0.5) + suction_mass_flow * (sf.H + (suction_velocity ** 2) * 0.5)

            print("Right hand", right_hand)
            print("Left hand" , left_hand)
            if math.fabs(right_hand - left_hand) <= 1000:
                return pf #Primary flow contains all relevant information
                break
            else:
                assumed_suction_pressure = assumed_suction_pressure - 100
                print('Assumed Pressure ', assumed_suction_pressure)            

ejector = Ejector()
flow_throat, velocity_throat = ejector.calculate_nozzle_flow()
print('Pressure in Nozzle Throat is ', flow_throat.P)

flow1, velocity_outlet = ejector.calculate_nozzle_outlet_flow(flow_throat, velocity_throat)
print('Calling first')

print('Pressure in Nozzle Outlet is ', flow1.P)
print('Mass flow rate is ', flow_throat.rho * ejector.nozzle_throat_area * velocity_throat)
mixing_flow = ejector.calculate_secondary_pressure(flow1, velocity_outlet)

Thanks so much for the quick reply! I'm not entirely sure what's going on, I've been examining the entropy values more closely, and they're negative.

For reference: 

Nozzle temp is 275, Nozzle pressure is 159000.0, Nozzle entropy is -389.13764658464726

Actual value should be 6.38 KJ/kgK, can you explain why I'm seeing these values? Am I not calculating entropy correctly?

CalebBell commented 6 years ago

Hi Rob, Entropy and enthalpy are always calculated with respect to a reference point. To compare values from different sources, one must always use only changes in values, not absolute values. REFPROP has much more accurate values of enthalpy than thermo; thermo has only a cubic EOS behind it, whereas theirs is a custom EOS just for ammonia. I'd hope it wouldn't be too much of a difference though.

The issue with the solver is indeed a phase change you are observing; the entropy is discontinuous at the phase change, as shown in this crude plot. image You seem comfortable with the library; I'd encourage you to use the HeatCapacityLiquid and HeatCapacityGas (use the T_dependent_property_integral and T_dependent_property_integral_over_T methods) along with the eos departure enthalpy/entropies, and build your own entropy/enthalpy reference frame you can fully understand and be confident in from the pieces. It's a great learning experience! Feel free to ask more questions. Sincerely, Caleb

CalebBell commented 3 years ago

This functionality should be quite reliable now with the new flash interface: https://thermo.readthedocs.io/thermo.flash.html#thermo.flash.FlashPureVLS A number of components are flashed with all sorts of different specifications to a very low tolerance. The implemented phases include cubic equations of state as well as ideal and activity coefficient based models: https://thermo.readthedocs.io/thermo.phases.html Hopefully no one else will need to go through your pain, Pheenixm! Sincerely, Caleb