more work on reactive distillation

benchmark/lactic_acid_purification_system.py

@@ -0,0 +1,144 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sat Mar 16 13:38:11 2024
+
+@author: cortespea
+"""
+import biosteam as bst
+from thermosteam.constants import R
+from math import exp
+try:
+ from .profile import register
+except:
+ def register(*args, **kwargs):
+ return lambda f: f
+
+__all__ = (
+ 'create_system_lactic_acid_purification',
+)
+
+@register(
+ 'lactic_acid_purification', 'Lactic acid purification',
+ 10, [2, 4, 6, 8, 10], 'LA\nsep.'
+)
+def create_system_lactic_acid_purification(alg='sequential modular'):
+ bst.settings.set_thermo(['Water', 'LacticAcid', 'ButylLactate', 'Butanol'], cache=True)
+
+ class Esterification(bst.KineticReaction):
+
+ def volume(self, stream): # kg of catalyst
+ rho_cat = 770 # kg / m3
+ liquid_volume = self.liquid_volume
+ catalyst_volume = 0.5 * liquid_volume
+ catalyst_mass = catalyst_volume * rho_cat
+ return catalyst_mass
+
+ def rate(self, stream): # kmol/kg-catalyst/hr
+ T = stream.T
+ # if T > 370: return 0 # Prevents multiple steady states
+ kf = 2.59e4 * exp(-5.340e4 / (R * T))
+ kr = 3.80e3 * exp(-5.224e4 / (R * T))
+ H2O, LA, BuLA, BuOH = stream.mol / stream.F_mol
+ return self.stoichiometry * 3600 * (kf * LA * BuOH - kr * BuLA * H2O) # kmol / kg-catalyst / hr
+
+ with bst.System(algorithm=alg) as sys:
+ feed = bst.Stream(
+ 'feed',
+ LacticAcid=4.174,
+ Water=5.470,
+ )
+ makeup_butanol = bst.Stream('makeup_butanol')
+ recycle_butanol = bst.Stream('recycle_butanol')
+ esterification_reflux = bst.Stream('esterification_reflux')
+ esterification = bst.MESHDistillation(
+ 'esterification',
+ ins=(feed, makeup_butanol, recycle_butanol, esterification_reflux), 
+ outs=('empty', 'bottoms', 'esterification_distillate'),
+ N_stages=17,
+ feed_stages=(1, 16, 16, 0),
+ stage_specifications={
+ 16: ('Boilup', 1),
+ 0: ('Boilup', 0),
+ },
+ liquid_side_draws={
+ 0: 1.0,
+ },
+ stage_reactions={
+ i: Esterification('LacticAcid + Butanol -> Water + ButylLactate', reactant='LacticAcid')
+ for i in range(1, 17)
+ },
+ maxiter=200,
+ LHK=('Butanol', 'ButylLactate'),
+ P=0.3 * 101325,
+ )
+ @esterification.add_specification(run=True)
+ def adjust_flow():
+ target = 5.85
+ makeup_butanol.imol['Butanol'] = max(target - recycle_butanol.imol['Butanol'], 0)
+
+ esterification_settler = bst.StageEquilibrium(
+ 'esterification_settler',
+ ins=(esterification-2), 
+ outs=(esterification_reflux, 'water_rich'),
+ phases=('L', 'l'),
+ top_chemical='Butanol',
+ )
+ water_distiller = bst.BinaryDistillation(
+ ins=esterification_settler-1, outs=('water_rich_azeotrope', 'water'),
+ x_bot=0.0001, y_top=0.2, k=1.2, Rmin=0.01,
+ LHK=('Butanol', 'Water'),
+ )
+ splitter = bst.Splitter(ins=water_distiller-1, split=0.5) # TODO: optimize split
+ hydrolysis_reflux = bst.Stream('hydrolysis_reflux')
+ hydrolysis = bst.MESHDistillation(
+ 'hydrolysis',
+ ins=(esterification-1, splitter-0, hydrolysis_reflux),
+ outs=('empty', 'lactic_acid', 'hydrolysis_distillate'),
+ N_stages=53,
+ feed_stages=(27, 50, 0),
+ stage_specifications={
+ 0: ('Boilup', 0),
+ 52: ('Boilup', 1),
+ },
+ liquid_side_draws={
+ 0: 1.0,
+ },
+ stage_reactions={
+ i: Esterification('LacticAcid + Butanol -> Water + ButylLactate', reactant='LacticAcid')
+ for i in range(1, 52) # It will run in reverse
+ },
+ P=101325,
+ LHK=('Butanol', 'LacticAcid'),
+ )
+
+ # @esterification.add_specification(run=True)
+ # def adjust_flow():
+ # target = 5.85
+ # makeup_butanol.imol['Butanol'] = max(target - recycle_butanol.imol['Butanol'], 0)
+
+ # Decanter
+ butanol_rich_azeotrope = bst.Stream('butanol_rich_azeotrope')
+ hydrolysis_settler = bst.StageEquilibrium(
+ 'settler',
+ ins=(hydrolysis-2, water_distiller-0, butanol_rich_azeotrope), 
+ outs=('butanol_rich_extract', hydrolysis_reflux),
+ phases=('L', 'l'),
+ top_chemical='Butanol',
+ T=310,
+ )
+
+ # Butanol purification
+ butanol_distiller = bst.BinaryDistillation(
+ ins=(hydrolysis_settler-0),
+ outs=(butanol_rich_azeotrope, recycle_butanol),
+ x_bot=0.0001, y_top=0.6, k=1.2, Rmin=0.01,
+ LHK=('Water', 'Butanol'),
+ )
+
+ return sys
+
+if __name__ == '__main__':
+ sys = create_system_lactic_acid_purification()
+ sys.flatten()
+ sys.diagram()
+ sys.simulate()

biosteam/_unit.py

@@ -1042,8 +1042,6 @@ def utility_cost(self) -> float:
  """Total utility cost [USD/hr]."""
  return self._utility_cost
 
-
-
  def mass_balance_error(self):
  """Return error in stoichiometric mass balance. If positive,
  mass is being created. If negative, mass is being destroyed."""

biosteam/units/distillation.py

@@ -2493,6 +2493,7 @@ def _init(self,
  method=None,
  inside_out=None,
  maxiter=None,
+ stage_specifications=None,
  ):
  if full_condenser: 
  if liquid_side_draws is None:
@@ -2501,10 +2502,11 @@ def _init(self,
  liquid_side_draws[0] = reflux / (1 + reflux)
  reflux = inf # Boil-up is 0
  self.LHK = LHK
- stage_specifications = {}
- stage_specifications[0] = ('Reflux', reflux)
- stage_specifications[-1] = ('Boilup', boilup)
-
+ if stage_specifications is None: stage_specifications = {}
+ if reflux is not None:
+ stage_specifications[0] = ('Reflux', reflux)
+ if boilup is not None:
+ stage_specifications[-1] = ('Boilup', boilup)
  super()._init(N_stages=N_stages, feed_stages=feed_stages,
  top_side_draws=vapor_side_draws, 
  bottom_side_draws=liquid_side_draws, 
@@ -2548,8 +2550,13 @@ def reflux(self, reflux):
 
  @property
  def boilup(self):
- name, value = self.stage_specifications[-1]
- if name == 'Boilup': return value
+ if -1 in self.stage_specifications:
+ name, value = self.stage_specifications[-1]
+ if name == 'Boilup': return value
+ else:
+ last_stage = self.N_stages - 1
+ name, value = self.stage_specifications[last_stage]
+ if name == 'Boilup': return value
  @boilup.setter
  def boilup(self, boilup):
  self.stage_specifications[-1] = ('Boilup', boilup)
@@ -2558,7 +2565,7 @@ def _setup(self):
  super()._setup()
  args = (self.N_stages, self.feed_stages, self.vapor_side_draws, 
  self.liquid_side_draws, self.use_cache, *self._ins, 
- self.partition_data, self.P,
+ self.partition_data, self.P, 
  self.reflux, self.boilup)
  if args != self._last_args:
  MultiStageEquilibrium._init(
@@ -2569,6 +2576,7 @@ def _setup(self):
  bottom_side_draws=self.liquid_side_draws,
  partition_data=self.partition_data, 
  stage_specifications=self.stage_specifications,
+ stage_reactions=self.stage_reactions,
  use_cache=self.use_cache, 
  )
  self._last_args = args
@@ -2668,18 +2676,24 @@ def update_liquid_holdup(self):
  for i in self.stage_reactions:
  partition = partitions[i]
  vapor, liquid = partition.outs
- rho_V = vapor.rho
- rho_L = liquid.rho
- active_area = area * (1 - partition.downcomer_area_fraction)
- Ua = vapor.get_total_flow('ft3/s') / active_area
- Ks = Ua * sqrt(rho_V / (rho_L - rho_V)) # Capacity parameter
+ if vapor.isempty():
+ Ks = 0
+ elif liquid.isempty():
+ partition.reaction.liquid_volume = 0
+ continue
+ else:
+ rho_V = vapor.rho
+ rho_L = liquid.rho
+ active_area = area * (1 - partition.downcomer_area_fraction)
+ Ua = vapor.get_total_flow('ft3/s') / active_area
+ Ks = Ua * sqrt(rho_V / (rho_L - rho_V)) # Capacity parameter 
  Phi_e = exp(-4.257 * Ks**0.91) # Effective relative froth density 
  Lw = 0.73 * diameter * 12 # Weir length [in] assuming Ad/A = 0.1
  # TODO: Compute weir length or other Ad/A
  qL = liquid.get_total_flow('gal/min')
  CL = 0.362 + 0.317 * exp(-3.5 * weir_height)
  hL = Phi_e * (hw + CL * (qL / (Lw * Phi_e)) ** b) # equivalent height of clear liquid holdup [in]
- partition.liquid_volume = hL * area * 0.00236155 # m3 
+ partition.reaction.liquid_volume = hL * area * 0.00236155 # m3 
 
  def estimate_diameter(self): # ft
  diameters = []