Recycled stream issue

[wenjung2]

Hi Yoel! Hope the codes below can reproduce the error.

When only simulating the MeOH_SynthesisReactor, it works well.

But when simulating system with more units and recycled streams (including complex heat exchange process as the attached figure), it gives:
*** RuntimeError: Failed to extrapolate gas heat capacity method 'HEOS_FIT' at T=nan K for component with CASRN '1333-74-0'.

I think some recycled stream caused this issue and maybe you have more quicker way to figure it out. Could you please help me check which stream causes this issue? Thanks!!

[_chemicals.py]

import thermosteam as tmo

#%% 
# Constants

_cal2joule=4.184

#%% 

chems=tmo.Chemicals([])

database_chemicals_dict={}
copied_chemicals_dict={}
defined_chemicals_dict={}

def chemical_database(ID, search_ID=None, phase=None, **kwargs):
    chemical = tmo.Chemical(ID,search_ID=search_ID, **kwargs)
    if phase:
        chemical.at_state(phase)
        chemical.phase_ref = phase
    chems.append(chemical)
    database_chemicals_dict[ID] = f'{ID}: {chemical.formula}/{chemical.MW}'
    return chemical

def chemical_copied(ID,ref_chemical,**data):
    chemical=ref_chemical.copy(ID)
    chems.append(chemical)
    for i,j in data.items():setattr(chemical,i,j)
    copied_chemicals_dict[ID]=f'{ID}:{chemical.formula}/{chemical.MW}'
    return chemical

def chemical_defined(ID,**kwargs):
    chemical=tmo.Chemical.blank(ID,**kwargs)
    chems.append(chemical)
    defined_chemicals_dict[ID]=f'{ID}:{chemical.formula}/{chemical.MW}'
    return chemical

# %% 

# =============================================================================
# Create chemical objects available in database
# =============================================================================

H2O = chemical_database('H2O')

O2 = chemical_database('O2', phase='g', Hf=0)
N2 = chemical_database('N2', phase='g', Hf=0)
H2 = chemical_database('H2', Hf=0)
CH4 = chemical_database('Methane')
CO = chemical_database('CarbonMonoxide', Hf=-26400*_cal2joule)
CO2 = chemical_database('CO2')

CH3OH = chemical_database('Methanol')

HCOOH = chemical_database('HCOOH')
C18H39N = chemical_database('C18H39N') # amine
C19H41NO2 = chemical_database('C19H41NO2')

Triphenylphosphine = chemical_database('Triphenylphosphine', search_ID='603-35-0', phase='s') # surrogate model for DCPE

DCPE = chemical_database('DCPE', search_ID='23743-26-2', phase='s') # FA catalyst; second component
for i in DCPE.get_missing_properties():
    if i == 'V':
        try:
            DCPE.copy_models_from(Triphenylphosphine, [i])
        except:
            pass

phase_change_chemicals = ['H2O', 'H2', 'CH4', 'CO', 'CO2', 'CH3OH', 'HCOOH'
                          'C18H39N', 'C19H41NO2']

for chem in chems:
    if chem.ID in phase_change_chemicals: pass
    elif chem.locked_state: pass
    else: 
        # Set phase_ref to avoid missing model errors
        if chem.phase_ref == 'g':
            chem.at_state('g')
        if chem.phase_ref == 'l':
            chem.at_state('l')

chems.compile()
tmo.settings.set_thermo(chems, cache=True)

chems.set_synonym('H2O', 'Water')
chems.set_synonym('CO2', 'CarbonDioxide')
chems.set_synonym('CO', 'CarbonMonoxide')
chems.set_synonym('Methanol', 'CH3OH')

[_units.py]

import numpy as np
import biosteam as bst
import thermosteam as tmo
from biosteam import Stream, Unit, BinaryDistillation
from biosteam.units import HXutility, Mixer, SolidsSeparator, Compressor
from biosteam.units.decorators import cost
from biosteam.units.design_tools import size_batch
from thermosteam import MultiStream
from biosteam.units.design_tools.geometry import cylinder_diameter_from_volume
from qsdsan._sanunit import SanUnit

Rxn = tmo.reaction.Reaction
ParallelRxn = tmo.reaction.ParallelReaction

#%% 
# =============================================================================
# H2 production through Liquid Alkaline (LA) electrolysis 
# =============================================================================
# Model H2 production by H2O => H2 + 0.5O2
class Electrolyzer(SanUnit): 
    _N_ins = 1
    _N_outs = 2
    
    def __init__(self, ID="", ins=None, outs=(),
                 stack_lifetime=10,
                 renewable_electricity_price=0.03, # currently possible in U.S. markets with plentiful  wind
                 **kwargs):
        bst.Unit.__init__(self, ID, ins, outs)
        self.stack_lifetime = stack_lifetime
        self.renewable_electricity_price = renewable_electricity_price
        
    def _run(self):
        water, = self.ins
        hydrogen, oxygen = self.outs
        hydrogen.phase = 'g'
        hydrogen.P = 30 * 101325 # Model H2 purification so that output is 99.99% purity and 30 bar
        oxygen.phase = 'g'
        hydrogen.imol['H2'] = water.F_mol
        oxygen.imol['O2'] = water.F_mol / 2
    
    def _cost(self): # based on https://www.osti.gov/biblio/2203367 P17
        self.power_utility.rate = self.outs[0].imass['H2'] * 122000 * 24 / 50000 # base case 50 MTD / 122 MW
        # self.power_utility.cost = self.renewable_electricity_price * self.power_utility.rate
        self.baseline_purchase_costs['Stack'] = self.power_utility.rate * 292 # process equipment, piping, valves, and  instrumentation including temperature, pressure, flow, and level indicators
        self.baseline_purchase_costs['Mechanical BOP'] = self.power_utility.rate * 280
        self.baseline_purchase_costs['Electrical BOP'] = self.power_utility.rate * 155 # for AC-> DC, transformer, power substation
        self.add_OPEX = self._calculate_replacement_cost()
    
    def _calculate_replacement_cost(self):
        stack_replacement_cost  = self.baseline_purchase_costs['Stack'] / self.stack_lifetime / (365*24) # convert from USD/year to USD/hour
        return stack_replacement_cost
    


    

# =============================================================================
# MeOH synthesis
# =============================================================================

# Reactor is modelled as an adiabatic ideal plug flow reactor (PFR)
class MeOH_SynthesisReactor(bst.units.design_tools.PressureVessel, bst.Unit):
    _N_ins = _N_outs = 1
    _N_heat_utilities = 0
    
    _units = {**bst.design_tools.PressureVessel._units,
              'Volume': 'm^3',
              'Catalyst weight': 'kg'}
    
    _F_BM_default = {**bst.design_tools.PressureVessel._F_BM_default,
                     'Catalyst': 1.0} # Cu/ZnO/Al2O3
    
    def __init__(self, ID="", ins=None, outs=(), thermo=None, *,
                 P=7600000,
                 tau=6/3600, # tau=5-7s to h; from 10.1016/j.jcou.2015.01.006
                 vessel_material='Stainless steel 304',
                 vessel_type='Vertical',
                 length_to_diameter=0.15/0.016, # 10.1016/j.jclepro.2013.06.008
                 catalyst_density=1775, # In kg/m^3; 10.1016/j.jclepro.2013.06.008
                 catalyst_price=13,
                 porosity=0.5,):
                 
            super().__init__(ID, ins, outs, thermo)
            self.P = P
            self.tau = tau
            self.WHSV = 1/self.tau
            self.vessel_material = vessel_material # Vessel material
            self.vessel_type = vessel_type # 'Horizontal' or 'Vertical'
            self.length_to_diameter = length_to_diameter #: Length to diameter ratio
            self.catalyst_density = catalyst_density
            self.catalyst_price = catalyst_price # In USD/kg
            self.porosity = porosity #  fraction of the bed volume occupied by void space
            self.reactions = ParallelRxn([
                # Reaction definition                 Reactant    Conversion
                Rxn('CO2 + 3H2 -> CH3OH + H2O',       'CO2',      0.210),
                Rxn('CO2 + H2 -> CO + H2O',           'CO2',      0.004),])
    
    def _default_vessel_type(self):
        return 'Vertical'
    
    def _run(self):
        feed, = self.ins
        assert feed.phase == 'g', 'feed phase must be a gas'
        effluent, = self.outs
        effluent.copy_like(feed)
        self.reactions.adiabatic_reaction(effluent)
        
    def _design(self):
        feed, = self.ins
        effluent, = self.outs
        length_to_diameter = self.length_to_diameter
        P = effluent.get_property('P', 'psi')
        results = self.design_results
        catalyst = feed.F_mass / self.WHSV
        results['Catalyst weight'] = catalyst
        results['Volume'] = volume = feed.F_vol * self.tau / self.porosity
        D = cylinder_diameter_from_volume(volume, length_to_diameter)
        L = length_to_diameter * D
        results.update(self._vessel_design(P, D, L))
        
    def _cost(self):
        D = self.design_results
        self.purchase_costs.update(
            self._vessel_purchase_cost(D['Weight'], D['Diameter'], D['Length'])
        )
        self.purchase_costs['Catalyst'] = self.catalyst_price * D['Catalyst weight']

[system.py]

import numpy as np
import thermosteam as tmo
import biosteam as bst
from biosteam import Flowsheet as F 
from biosteam import main_flowsheet
from biosteam import units, Stream, SystemFactory
from biosteam.process_tools import UnitGroup
from biorefineries.ccu._chemicals import chems
from biorefineries.ccu import _units

F = bst.Flowsheet('CO2_methane')
bst.main_flowsheet.set_flowsheet(F)

tmo.settings.set_thermo(chems, cache=True)

#%%
# =============================================================================
# 1. MeOH production; based on 
# Design and simulation of a methanol production plant from CO2 hydrogenation
# =============================================================================
#add catalyst
water_stream_1 = Stream('water_stream_1',
                      Water=1,
                      phase='l',
                      total_flow=10000,
                      units='kg/hr')

CO2_stream_1 = Stream('CO2_stream_1',
                    CO2=1,
                    phase='g',
                    total_flow=88000,
                    units='kg/hr')

# h2 = Stream('h2',
#             H2=1,
#             phase='g',
#             T=25+273.15,
#             P=30*101325,
#             total_flow=12100,
#             units='kg/hr')
                    

# CO2 compressing
ks = [bst.units.IsentropicCompressor(P=3*101325), 
      bst.units.IsentropicCompressor(P=8*101325),
      bst.units.IsentropicCompressor(P=26.3*101325)]

hxs = [bst.units.HXutility(T=T) for T in [38+273.15, 38+273.15, 38+273.15]]

C1101 = units.MultistageCompressor('C1101', ins=CO2_stream_1, outs='', compressors=ks, hxs=hxs)

C1102 = units.IsentropicCompressor('C1102', ins=C1101-0, outs='', P=78*101325)

# Water electrolysis
R1101 = _units.Electrolyzer('R1101', ins=water_stream_1, outs=('hydrogen','oxygen'))
@R1101.add_specification(run=True)
def adjust_water_flow():
    C1101.run()
    R1101.ins[0].F_mol = C1101.ins[0].F_mol * 3 # H2:CO2 = 3:1

# H2 compressing to same pressure
ks_2 = [bst.units.IsentropicCompressor(P=78*101325)]

hxs_2 = [bst.units.HXutility(T=164+273.15)]

C1103 = units.MultistageCompressor('C1103', ins=R1101-0, outs='', compressors=ks_2, hxs=hxs_2)

M1101 = units.Mixer('M1101', ins=(C1102-0, C1103-0), outs='', rigorous=True)

M1102 = units.Mixer('M1102', ins=(M1101-0, ''), outs='')

H1101 = units.HXprocess('H1101', ins=(M1102-0, ''), outs=('', ''))

# H1101_1 = units.HXutility('H1101_1', ins=H1101-0, outs='', T=210+273.15)

R1102 = _units.MeOH_SynthesisReactor('R1102', ins=H1101-0, outs='product')

S1101 = units.Splitter('S1101', ins=R1102-0, outs=(1-H1101, ''), split=0.65)

V1101 = units.IsenthalpicValve('V1101', ins=H1101-1, outs='', P=73.6*101325)

M1103 = units.Mixer('M1103', ins=(V1101-0, ''), outs='')

H1102 = units.HXutility('H1102', ins=M1103-0, outs='', T=35+273.15, rigorous=True)

S1102 = units.PhaseSplitter('S1102', ins=H1102-0, outs=('gas', 'condensed_water_and_methanol'))

S1103 = units.Splitter('S1103', ins=S1102-0, outs=('purge', 'recycled'), split=0.01)

C1104 = units.IsentropicCompressor('C1104', ins=S1103-1, outs=1-M1102, P=78*101325)

# For condensed water and methanol in S1102 (composed of methanol, water and residual dissolved gases)
V1102 = units.IsenthalpicValve('V1102', ins=S1102-1, outs='', P=10*101325)

V1103 = units.IsenthalpicValve('V1103', ins=V1102-0, outs='', P=1.2*101325)

F1101 = units.SplitFlash('F1101', ins=V1103-0, outs=('gas', ''), T=22+273.15, P=1.2*101325,
                         split=dict(CO2=0.999, H2=0.999))
                                                                                                       

H1103 = units.HXprocess('H1103', ins=(S1101-1, F1101-1), outs=(1-M1103, ''))

D1101 = units.BinaryDistillation('D1101', ins=H1103-1, outs=('gas_MEOH', 'bottom_water'),
                                 LHK=('Methanol', 'Water'),
                                 Lr=0.9999, Hr=0.9999, k=2,
                                 is_divided=True)


sys = bst.main_flowsheet.create_system('MeOH_sys')
sys.set_tolerance(method='fixedpoint')
sys.maxiter = 600 # 200 will not lead to convergence
sys.empty_recycles()
sys.simulate()
sys.diagram()

[yoelcortes]

@wenjung2, thanks for posting the code. Could you also simplify it to a minimal reproducible example?

Since the problem is just the MeOH synthesis reactor, we just need the MeOH systhesis code and the feed stream. You can run the code that creates the error, then copy the inlet feed stream values. The other code is not relevant.

Considering the test below, it should work if the stream conditions are the same. You may want to add a preheater bst.HXutility(..., T=210 + 273.15, V=1, rigorous=False) to force the same thermodynamic conditions.

import biosteam as bst
bst.settings.set_thermo(['CO2', 'H2', 'CH3OH', 'H2O'])
rxn = bst.Rxn('CO2 + 3H2 -> CH3OH + H2O', reactant='CO2')
s = bst.Stream(CO2=1, H2=3, phase='g', P=7.8e6, T=210 + 273.15)
s.show()
rxn.adiabatic_reaction(s)
s.show()

Output:

Stream: s1
phase: 'g', T: 483.15 K, P: 7.8e+06 Pa
flow (kmol/hr): CO2  1
                H2   3
Stream: s1
phase: 'g', T: 990.25 K, P: 7.8e+06 Pa
flow (kmol/hr): CH3OH  1
                H2O    1

Thanks,

[wenjung2]

Thanks Yoel! Currently I am fully occupied by other work but I will look into it later.