fix infeasible default aerobic bioreactor values; add tests

Author: yoelcortes

biosteam/__init__.py

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

biosteam/units/abstract_stirred_tank_reactor.py

@@ -51,6 +51,8 @@ class AbstractStirredTankReactor(PressureVessel, Unit, isabstract=True):
         Fraction of working volume over total volume. Defaults to 0.8.
     V_max :
         Maximum volume of a reactor [m3]. Defaults to 355.
+    length_to_diameter :
+        Length to diameter ratio of bioreactor.
     kW_per_m3 : 
         Power usage of agitator. Defaults to 0.985 [kW / m3] converted from 
         5 hp/1000 gal as in [1]_, for liquid–liquid reaction or extraction.

biosteam/units/aerated_bioreactor.py

@@ -75,6 +75,8 @@ class AeratedBioreactor(AbstractStirredTankReactor):
         Operating pressure [Pa].
     V_wf : 
         Fraction of working volume over total volume. Defaults to 0.8.
+    length_to_diameter :
+        Length to diameter ratio of bioreactor.
     V_max :
         Maximum volume of a reactor [m3]. Defaults to 355.
     kW_per_m3 : 
@@ -205,15 +207,23 @@ def _init(
         AbstractStirredTankReactor._init(self, **kwargs)
         self.theta_O2 = theta_O2 # Average concentration of O2 in the liquid as a fraction of saturation.
         self.Q_O2_consumption = Q_O2_consumption # Forced duty per O2 consummed [kJ/kmol].
-        self.optimize_power = True if optimize_power is None else optimize_power
         if design is None: 
-            design = 'Stirred tank'
+            if self.kW_per_m3 == 0:
+                design = 'Bubble column'
+            else:
+                design = 'Stirred tank'
         elif design not in aeration.kLa_method_names:
             raise ValueError(
                 f"{design!r} is not a valid design; only "
                 f"{list(aeration.kLa_method_names)} are valid"
             )
         self.design = design
+        self.optimize_power = (
+            design == 'Stirred tank' 
+            if optimize_power is None else
+            optimize_power
+        )
+        self.design = design
         if method is None:
             method = self.default_methods[design]
         if (key:=(design, method)) in aeration.kLa_methods:
@@ -271,7 +281,10 @@ def get_agitation_power(self, kLa):
             self.superficial_gas_flow = U = F / A # m / s
             return aeration.P_at_kLa_Riet(kLa, V, U, **self.kLa_kwargs)
         else:
-            raise NotImplementedError('kLa method has not been implemented in BioSTEAM yet')
+            raise NotImplementedError(
+                'cannot solve for the required agitation power using the '
+               f'kLa method {self.kLa!r} in BioSTEAM yet'
+            )
 
     def _get_duty(self):
         if self.Q_O2_consumption is None:
@@ -368,11 +381,31 @@ def air_flow_rate_objective(O2):
                 return OUR - self.get_OTR()
 
             f = air_flow_rate_objective
+            x0 = OUR
             y0 = air_flow_rate_objective(OUR)
-            if y0 <= 0.: # Correlation is not perfect and special cases lead to OTR > OUR
-                return
-            flx.IQ_interpolation(f, x0=OUR, x1=10 * OUR, 
-                                 y0=y0, ytol=1e-3, xtol=1e-3)
+            
+            # Correlation is not perfect and special cases lead to OTR > OUR
+            if y0 <= 0.: return 
+            
+            f = air_flow_rate_objective
+            x1 = 10 * OUR
+            y1 = air_flow_rate_objective(x1)
+            
+            # It is possible an infinite flow rate of air is not enough to 
+            # satisfy mass transfer if the titer is super high or gas O2 concentration
+            # is too low.
+            if y1 > 0: 
+                raise RuntimeError(
+                    'bioreactor conversion/titer cannot be satisfied; '
+                    'even an infinite flow rate of gas is not enough'
+                )
+            
+            # There is a known solution because y1 < 0 < y0
+            flx.IQ_interpolation(
+                f, x0=x0, x1=x1, 
+                y0=y0, y1=y1, 
+                ytol=1e-3, xtol=1e-3
+            )
 
     def _run_reactions(self, effluent):
         self.reactions.force_reaction(effluent)
@@ -422,7 +455,7 @@ def get_OTR(self):
         return OTR
 
     def _inlet_air_pressure(self):
-        AbstractStirredTankReactor._design(self)
+        AbstractStirredTankReactor._design(self, size_only=True)
         liquid = bst.Stream(None, thermo=self.thermo)
         liquid.mix_from([i for i in self.ins if i.phase != 'g'], energy_balance=False)
         liquid.copy_thermal_condition(self.outs[0])
@@ -433,14 +466,7 @@ def _inlet_air_pressure(self):
     def _design(self):
         AbstractStirredTankReactor._design(self)
         if self.air.isempty(): return
-        liquid = bst.Stream(None, thermo=self.thermo)
-        liquid.mix_from([i for i in self.ins if i.phase != 'g'], energy_balance=False)
-        liquid.copy_thermal_condition(self.outs[0])
-        rho = liquid.rho
-        length = self.get_design_result('Length', 'm') * self.V_wf
-        compressor = self.compressor
-        compressor.P = g * rho * length + 101325
-        compressor.simulate()
+        self.compressor.simulate()
         air_cooler = self.air_cooler
         air_cooler.T = self.T
         air_cooler.simulate()
@@ -508,6 +534,8 @@ class GasFedBioreactor(AbstractStirredTankReactor):
         Operating pressure [Pa].
     V_wf : 
         Fraction of working volume over total volume. Defaults to 0.8.
+    length_to_diameter :
+        Length to diameter ratio of bioreactor.
     V_max :
         Maximum volume of a reactor [m3]. Defaults to 355.
     kW_per_m3 : 
@@ -762,7 +790,8 @@ def _run(self):
         vent.T = effluent.T = self.T
         vent.empty()
         vent.phase = 'g'
-        if not self.titer:
+        if not self.titer: 
+            # Titer is given by the mass transfer.
             liquid_feeds = [i for i in self.ins if i.phase != 'g']
             T = self.T
             P = self._inlet_gas_pressure()
@@ -776,6 +805,7 @@ def _run(self):
             vent.empty()
             self._run_vent(vent, effluent)
         elif variable_gas_feeds:
+            # Solve gas flow rates to meet titer.
             liquid_feeds = [i for i in self.ins if i.phase != 'g']
             effluent.mix_from(liquid_feeds, energy_balance=False)
             T = self.T
@@ -806,7 +836,13 @@ def load_flow_rates(F_feeds):
 
             baseline_feed = bst.Stream.sum(self.normal_gas_feeds, energy_balance=False)
             baseline_flows = baseline_feed.get_flow('mol/s', self.gas_substrates)
-            bounds = np.array([[max(1.01 * SURs[i] - baseline_flows[i], 1e-6), 10 * SURs[i]] for i in index])
+            
+            # Bounds must meet substrate uptake rate (minimally).
+            # At most, 10x the substrate uptake rate as an abritrarily high number.
+            bounds = np.array([
+                [max(1.01 * SURs[i] - baseline_flows[i], 1e-6), 10 * SURs[i]] 
+                for i in index
+            ])
             if self.optimize_power:
                 def total_power_at_substrate_flow(F_substrates):
                     load_flow_rates(F_substrates)
@@ -837,8 +873,13 @@ def gas_flow_rate_objective(F_substrates):
                     bounds = bounds.T
                     results = least_squares(f, 1.2 * SURs, bounds=bounds, ftol=SURs.min() * 1e-6)
                 self._results = results
+                if not results.success:
+                    raise RuntimeError(
+                        'bioreactor conversion/titer could not be satisfied'
+                    )
                 load_flow_rates(results.x / x_substrates)
-        else:
+        else: 
+            # Titer given, must adjust liquid flows to meet mass transfer.
             try:
                 feed, = [i for i in self.ins if i.phase != 'g']
             except:
@@ -871,8 +912,6 @@ def f(vent_flow_rates):
                 return effluent.imass[product] / effluent.ivol['Water'] - titer
 
             flx.IQ_interpolation(liquid_flow_rate_objective, 0.05 * F_liquid_max, F_liquid_max, ytol=1e-3)
-        # self.show()
-        # breakpoint()
 
     def _run_reactions(self, effluent):
         data = effluent.get_data()

biosteam/units/anaerobic_bioreactor.py

@@ -41,6 +41,8 @@ class AnaerobicBioreactor(AbstractStirredTankReactor):
         Operating pressure [Pa].
     V_wf : 
         Fraction of working volume over total volume. Defaults to 0.8.
+    length_to_diameter :
+        Length to diameter ratio of bioreactor.
     V_max :
         Maximum volume of a reactor [m3]. Defaults to 355.
     kW_per_m3 : 

biosteam/units/single_phase_reactor.py

@@ -43,6 +43,8 @@ class SinglePhaseReactor(AbstractStirredTankReactor):
         Operating pressure [Pa].
     V_wf : 
         Fraction of working volume over total volume. Defaults to 0.8.
+    length_to_diameter :
+        Length to diameter ratio of bioreactor.
     V_max :
         Maximum volume of a reactor [m3]. Defaults to 355.
     kW_per_m3 : 

setup.py

@@ -11,7 +11,7 @@
     name='biosteam',
     packages=['biosteam'],
     license='MIT',
-    version='2.53.3',
+    version='2.53.4',
     description='The Biorefinery Simulation and Techno-Economic Analysis Modules',
     long_description=open('README.rst', encoding='utf-8').read(),
     author='Yoel Cortes-Pena',

tests/test_aerobic_bioreactor.py

@@ -0,0 +1,136 @@
+# -*- coding: utf-8 -*-
+# BioSTEAM: The Biorefinery Simulation and Techno-Economic Analysis Modules
+# Copyright (C) 2020-2023, Yoel Cortes-Pena <yoelcortes@gmail.com>
+# 
+# This module is under the UIUC open-source license. See 
+# github.com/BioSTEAMDevelopmentGroup/biosteam/blob/master/LICENSE.txt
+# for license details.
+"""
+"""
+import biosteam as bst
+from biorefineries.sugarcane import chemicals
+from numpy.testing import assert_allclose
+import pytest
+
+def test_jacketed_bubble_column():
+    # Get arguments ready
+    bst.settings.set_thermo(chemicals)
+    feed = bst.Stream(
+        Water=1.20e+05,
+        Glucose=2.5e+04,
+        units='kg/hr',
+        T=32+273.15
+    )
+    rxn = bst.Rxn(
+        'Glucose + O2 -> H2O + CO2', reactant='Glucose', X=0.5,
+        correct_atomic_balance=True
+    ) 
+    kwargs = dict(
+        ins=[feed, bst.Stream('air', phase='g')],
+        outs=('vent', 'product'), tau=12, V_max=500,
+        design='Bubble column', reactions=rxn,
+    )
+    
+    # Must raise error if optimized power because default 
+    # kLa method does not depend on agitation.
+    R1 = bst.AeratedBioreactor(optimize_power=True, **kwargs)
+    with pytest.raises(NotImplementedError):
+        R1.simulate()
+
+    R1 = bst.AeratedBioreactor(**kwargs)
+    assert not R1.optimize_power, 'optimize power should default to false'
+    
+    # Test design results
+    R1.simulate()
+    expected_results = {
+        'Reactor volume': 453.1780985777661, 
+        'Batch time': 18.0, 
+        'Loading time': 2.9999999999999996, 
+        'Residence time': 12, 
+        'Length': 56.81446196104029, 
+        'Diameter': 18.93815398701343, 
+        'Weight': 73673.55, 
+        'Wall thickness': 0.39277612951101165
+    }
+    for name, expected in expected_results.items():
+        assert_allclose(R1.design_results[name], expected)
+
+    assert R1.design_results['Vessel type'] == 'Vertical', 'vessel type must be vertial'
+
+    # Test gas flow rate sensitivity to aspect ratio
+    air_flow_rates = []
+    for L2D in (1, 3, 5):
+        R1.length_to_diameter = L2D
+        R1.simulate()
+        air_flow_rates.append(R1.air.F_mass)
+    
+    assert_allclose(
+        air_flow_rates, 
+        [378670.3313278283, 
+         188239.2629457292, 
+         137890.44337525323]
+    )
+    
+    # Must raise error if titer/yield is not feasible (mass transfer limitation)
+    R1.length_to_diameter = 0.01
+    with pytest.raises(RuntimeError): R1.simulate()
+
+
+def test_jacketed_stirred_tank():
+    # Get arguments ready
+    bst.settings.set_thermo(chemicals)
+    feed = bst.Stream(
+        Water=1.20e+05,
+        Glucose=2.5e+04,
+        units='kg/hr',
+        T=32+273.15
+    )
+    rxn = bst.Rxn(
+        'Glucose + O2 -> H2O + CO2', reactant='Glucose', X=0.5,
+        correct_atomic_balance=True
+    ) 
+    kwargs = dict(
+        ins=[feed, bst.Stream('air', phase='g')],
+        outs=('vent', 'product'), tau=12, V_max=500,
+        reactions=rxn,
+    )
+    R1 = bst.AeratedBioreactor(**kwargs)
+    assert R1.design == 'Stirred tank'
+    assert R1.optimize_power, 'optimize power should default to True'
+    
+    # Test design results
+    R1.simulate()
+    expected_results = {
+        'Reactor volume': 453.1780985777661, 
+        'Batch time': 18.0, 
+        'Loading time': 2.9999999999999996, 
+        'Residence time': 12, 
+        'Length': 56.81446196104029, 
+        'Diameter': 18.93815398701343, 
+        'Weight': 73673.55, 
+        'Wall thickness': 0.39277612951101165
+    }
+    for name, expected in expected_results.items():
+        assert_allclose(R1.design_results[name], expected)
+
+    assert R1.design_results['Vessel type'] == 'Vertical', 'vessel type must be vertial'
+
+    # Test gas flow rate sensitivity to aspect ratio
+    air_flow_rates = []
+    for L2D in (1, 3, 5):
+        R1.length_to_diameter = L2D
+        R1.simulate()
+        air_flow_rates.append(R1.air.F_mass)
+    
+    assert_allclose(
+        air_flow_rates, 
+        [222037.13645703293, 
+         140467.1502394925, 
+         116058.21484691999]
+    )
+    
+
+if __name__ == '__main__':
+    test_jacketed_bubble_column()
+    test_jacketed_stirred_tank()
+