Recycled stream issue #185

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')

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.201501006
 vessel_material='Stainless steel 304',
 vessel_type='Vertical',
 length_to_diameter=0.15/0.016, # 10.1016/j.jclepro.201306008
 catalyst_density=1775, # In kg/m^3; 10.1016/j.jclepro.201306008
 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']

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()