0.52.16 release; add tutorial on distillation

Author: yoelcortes

biosteam/__init__.py

@@ -13,7 +13,7 @@
 
 """
 from __future__ import annotations
-__version__ = '2.52.15'
+__version__ = '2.52.16'
 
 #: Chemical engineering plant cost index (defaults to 567.5 at 2017).
 CE: float = 567.5 

biosteam/units/distillation.py

@@ -431,6 +431,55 @@ def _load_components(self, partial_condenser, condenser_thermo, reboiler_thermo)
             self.reboiler-0, ('boilup', self-1), thermo=reboiler_thermo,
         )
 
+    def to_rigorous_column(self, dP=None, kwargs_only=False, print_code=False, **kwargs):
+        design = self.design_results
+        feed_stage = int(design['Theoretical feed stage'])
+        N_stages = int(design['Theoretical stages'])
+        R = design['Reflux']
+        flow_split = self.outs[1].F_mol / self.ins[0].F_mol
+        if dP is None:
+            dP = 690.0 # Heuristically, pressure drops are between 0.35 to 1.03 kPa
+        if dP:
+            P_reboiler = self.outs[1].P
+            P = [P_reboiler + i * dP for i in range(N_stages)]
+        else:
+            P = self.P
+        kwargs = dict(
+            N_stages=N_stages, 
+            feed_stages=[feed_stage],
+            LHK=self.LHK,
+            stage_specifications={
+                0: ('Reflux', R), 
+                -1: ('Flow', flow_split),
+            },
+            P=P,
+            **kwargs,
+        )
+        if print_code:
+            kwargs_code = kwargs.copy()
+            kwargs_code['P'] = f'[{self.P} + i*{dP} for i in range({N_stages})]'
+            kwargs_code['stage_specifications'] = {
+                0: ('Reflux', round(R, 3)), 
+                -1: ('Flow', round(flow_split, 5)),
+            }
+            code = (
+                f'MESHDistillation(\n'
+                f"    ins=[{', '.join([i.ID for i in self.ins])}],\n"
+                f"    outs=[{', '.join([i.ID for i in self.outs])}]"
+            )
+            code += bst.repr_kwargs(kwargs_code, dlim=',\n    ')
+            code += '\n)'
+            print(code)
+        if kwargs_only: 
+            return kwargs
+        return MESHDistillation(
+            None,
+            ins=self.ins, 
+            outs=self.outs,
+            thermo=self.thermo,
+            **kwargs
+        )
+    
     @property
     def distillate(self):
         return self.outs[0]
@@ -866,36 +915,33 @@ def _run_condenser_and_reboiler(self):
         liq = reboiler.ins[0]
         liq.mix_from([bottoms_product, boilup], energy_balance=False)
         liq.phase = 'l'
-        liq_T = liq.bubble_point_at_P().T
-        if liq_T > bp.T: liq_T = bp.T - 0.1
-        liq.T = liq_T
-        self.pump.ins[0].copy_like(liq)
-        self.pump.simulate()
+        liq.T = bp.T
         self.bottoms_split.simulate()
         if self._partial_condenser: self.reflux_drum.simulate()
 
     def _simulate_components(self): 
         reboiler = self.reboiler
         condenser = self.condenser
         Q_condenser = condenser.outs[0].H - condenser.ins[0].H
+        condenser_kwargs = dict(duty=Q_condenser)
         H_out = self.H_out
         H_in = self.H_in
         Q_overall_boiler =  H_out - H_in - Q_condenser
-        Q_boiler = reboiler.outs[0].H - reboiler.ins[0].H
-        if Q_boiler < Q_overall_boiler:
-            liquid = reboiler.ins[0]
-            H_out_boiler = reboiler.outs[0].H
-            try:
-                liquid.H = H_out_boiler - Q_overall_boiler
-            except:
-                liquid.phase = 'l'
-                boiler_kwargs = dict(duty=Q_boiler)                
-            else:
-                boiler_kwargs = dict(duty=Q_overall_boiler)
-            condenser_kwargs = dict(duty=Q_condenser)
+        H_out_boiler = reboiler.outs[0].H
+        liquid_in = reboiler.ins[0]
+        vapor = reboiler.outs[0]
+        try:
+            liquid_in.H = H_out_boiler - Q_overall_boiler
+        except:
+            liq_T = liquid_in.bubble_point_at_P().T
+            if liq_T > vapor.T: liq_T = vapor.T - 0.1
+            vapor.T = liq_T
+            Q_boiler = reboiler.outs[0].H - reboiler.ins[0].H
+            boiler_kwargs = dict(duty=Q_boiler)     
         else:
-            boiler_kwargs = dict(duty=Q_boiler)
-            condenser_kwargs = dict(duty=Q_condenser)
+            boiler_kwargs = dict(duty=Q_overall_boiler)
+        self.pump.ins[0].copy_like(liquid_in)
+        self.pump.simulate()
         reboiler.simulate(
             run=False,
             design_kwargs=boiler_kwargs,

biosteam/units/stage.py

@@ -2071,12 +2071,14 @@ class MultiStageEquilibrium(Unit):
     damping = 0 # Damping factor; defined as x_(i+1) = damping * x_i + (1 - damping) * f(x_i)
     minimum_residual_reduction = 0.25 # Minimum fractional reduction in residual for simulation.
     iteration_memory = 10 # Length of recorded iterations.
+    preconditioning_tolerance = 1e-3
+    preconditioning_relative_tolerance = 1e-3
     inside_maxiter = 100
     default_max_attempts = 5
     default_maxiter = 100
-    default_optimize_result = True
-    default_tolerance = 1e-5
-    default_relative_tolerance = 1e-4
+    default_optimize_result = False
+    default_tolerance = 1e-6
+    default_relative_tolerance = 1e-6
     default_algorithms = ('inside out', 'phenomena', 'simultaneous correction', 'sequential modular')
     default_inside_loop_algorithm = 'simultaneous correction'
     decomposition_algorithms = {
@@ -2258,7 +2260,7 @@ def _init(self,
             N_given = len(outs)
             N_outs = 2 + N_side_draws 
             top_mark = 2 + len(top_side_draws)
-            if stage_specifications[0] == ('Reflux', inf) and N_given < N_outs:
+            if stage_specifications.get(0) == ('Reflux', inf) and N_given < N_outs:
                 # Empty vapor is not given (for full condenser scenario)
                 tsd_iter = iter(outs[2:top_mark])
                 bsd_iter = iter([outs[0], *outs[top_mark:]])
@@ -2386,6 +2388,9 @@ def _init(self,
 
     # %% Optimization
 
+    def residual(self):
+        return self._objective(self._get_point())
+    
     def _get_point(self, x=None):
         if x is None: x = np.zeros(self._point_shape)
         for i, stage in enumerate(self.stages): 
@@ -2520,12 +2525,12 @@ def _run(self):
         x = self.hot_start()
         algorithms = self.algorithms
         methods = self.methods
-        optimize_result = self.optimize_result
+        optimize_result = self.optimize_result and 'K' not in self.partition_data
         f = self._iter
         analysis_mode = self._convergence_analysis_mode
         maxiter = self.maxiter
-        xtol = 100 * self.tolerance if optimize_result else self.tolerance
-        rtol = 100 * self.relative_tolerance if optimize_result else self.relative_tolerance
+        xtol = self.preconditioning_tolerance if optimize_result else self.tolerance
+        rtol = self.preconditioning_relative_tolerance if optimize_result else self.relative_tolerance
         self.iter = 0
         for n in range(self.max_attempts):
             self.attempt = n
@@ -2559,7 +2564,7 @@ def _run(self):
                     break
             else: continue
             break
-        if optimize_result and 'K' not in self.partition_data: x = self._simultaneous_correction(x, 'hybr')
+        if optimize_result: x = self._simultaneous_correction(x, 'hybr')
         self._set_point(x)
         # Last simulation to force mass balance
         self.update_mass_balance()

docs/tutorial/Distillation.ipynb

@@ -11,12 +11,12 @@
     name='biosteam',
     packages=['biosteam'],
     license='MIT',
-    version='2.52.15',
+    version='2.52.16',
     description='The Biorefinery Simulation and Techno-Economic Analysis Modules',
     long_description=open('README.rst', encoding='utf-8').read(),
     author='Yoel Cortes-Pena',
     install_requires=['IPython>=7.9.0',
-                      'thermosteam>=0.52.12', 
+                      'thermosteam>=0.52.16', 
                       'graphviz>=0.17',
                       'numpoly==1.2.13',
                       'chaospy==4.3.15',

setup.py

@@ -293,7 +293,7 @@ def test_distillation():
         N_stages=10,
         ins=[hot_extract],
         feed_stages=[5],
-        outs=['', 'bottoms_product', 'distillate'],
+        outs=['distillate', 'bottoms_product'],
         full_condenser=True,
         reflux=1.0,
         boilup=3.5,
@@ -302,9 +302,8 @@ def test_distillation():
     )
     distillation.simulate()
     flows = [
-        [0.0, 0.0, 0.0],
+        [22.055922550575424, 0.1438997998268773, 87.15249997005776],
         [0.11207744942457795, 4.350500200173122, 22.35850002994225],
-        [22.055922550575424, 0.1438997998268773, 87.15249997005776]
     ]
     for i, j in zip(distillation.outs, flows):    
         assert_allclose(i.mol, j, rtol=1e-3, atol=1e-3)