SanUnit (advanced)¶
Click the badge below to try this tutorial interactively in your browser:
You can also run this tutorial inGoogle Colab. It takes a one-time setup per session: follow theColab instructions.
Prepared by:
Learning objectives. After this tutorial, you will be able to:
Subclass
SanUnitto create a custom unit operationImplement
_run,_design, and_costmethodsIntegrate the new unit into a
System
Prerequisites: 4. SanUnit (basic)
Covered topics:
Basic structure of SanUnit subclasses
Making a simple AerobicReactor
Other convenient features
Companion video. A walkthrough of this tutorial is available on YouTube, presented by Hannah Lohman. Recorded against
QSDsanv1.2.0. The concepts still apply, but if the code on screen differs from this notebook, follow the notebook.
Setup¶
Import QSDsan and confirm the installed version.
[1]:
import qsdsan as qs
print(f'This tutorial was made with qsdsan v{qs.__version__}.')
This tutorial was made with qsdsan v1.5.3.
1. Basic structure of SanUnit subclasses¶
Building a custom unit means subclassing SanUnit. Assuming you are familiar with the topics covered in the previous tutorial on SanUnit, we can now learn the specifics of creating subclasses.
New to Python classes (the building blocks for subclassing)? Expand the aside below for a refresher, or see the resources in the tutorials index.
Python Aside: classes, methods, attributes, and properties (click to expand)
Classes and instances. A class bundles data and functions; you create instances from it.
class Apple:
kind = 'fruit' # class attribute (shared by all instances)
def __init__(self, name, color):
self.name = name # instance attributes (per object)
self.color = color
def introduce(self): # instance method (takes self)
print(f'{self.name} is {self.color}')
gala = Apple('Gala', 'red')
gala.introduce() # Gala is red
Method kinds. Instance methods take self; a @classmethod takes cls (e.g. Components.load_default); a @staticmethod takes neither. The @ is a decorator — a wrapper that adds behavior.
Properties. A property looks like an attribute but is computed by getter/setter functions, so you can validate or protect a value (omit the setter to make it read-only):
class Reactor:
@property
def conversion(self): # getter
return self._conversion
@conversion.setter
def conversion(self, i): # setter
if not 0 <= i <= 1:
raise AttributeError('conversion must be in [0, 1]')
self._conversion = i
Subclassing reuses a parent’s behavior: class Apple2(Apple): ... inherits introduce. Subclassing SanUnit is exactly how you build a custom unit — which is what the rest of this tutorial does.
1.1. Fundamental methods¶
SanUnit itself is a subclass of BioSTEAM’s Unit class. For more details, check out BioSTEAM’s documentation.
In addition to the __init__ method for initialization, all SanUnit objects have three most fundamental methods (they all start with _, as users typically don’t interact with them):
_run, which is used for mass and energy calculation within a unit operation (e.g., if you have an anaerobic reactor that will convert 80% of the organics in the influent, you’ll want to put it in_runThere is also a
runmethod that will call the_runmethod and anyspecificationfunctions you define, but we will skip it for now
_design, this method will be called after_run(when you have aSystemwith multiple units, then_designwill only be called after all units within the system have converged). The_designmethod contains algorithms that are used in designing the unit (e.g., volume, dimensions, materials)Material inventories calculated in
_designare usually stored in thedesign_resultsdict of the unit with keys being names (str) of the item and values being quantities of the itemAll entries in the
design_resultsdict should have corresponding entries in the_unitsdict to provide units of measure to the values in thedesign_resultsdict
_cost, which will be called after_designto get the cost of the unit, it may leverage the inventories calculated in_runor_designThe baseline purchase cost of each inventory item is usually stored in the
baseline_purchase_costsdict, and installed cost of this item will be calculated using different factors (F_BM,F_D,F_P, andF_Mfor bare module, design, pressure, and material factors, they are all dict).Each factor defaults to 1 if not given, but if a cost item is missing from
F_BM, BioSTEAM still costs it (bare-module factor 1) while emitting aRuntimeWarning. Declare the factor to silence the warning and be explicit: per item with a class-level_F_BM_defaultdict (e.g.,_F_BM_default = {'Tank': 2}) or by settingself.F_BM['Tank'] = 2, or set every item to 1 at once with theF_BM_default=1initialization argument (handy when a unit’s purchase costs already equal its installed costs).
These methods will be aggregated into a simulate function that will call these methods (and do some other minor stuff).
Note: You do NOT need to use all of these methods, and you do not need to strictly follow the functionalities above. For example, you can put cost algorithms in _design or even _run, but the latter will be strongly discouraged unless you have a good reason, as the cost algorithms will be run a lot of times when the System is trying to converge, which adds unnecessary overheads.
Heads up: a unit with no ``_run``. If a unit has exactly one inlet and one outlet and you do not define _run, It is treated as a static pass-through and links the outlet to the inlet, so the two share the same flow data (re-simulating can then drain the feed). That is handy for a genuine pass-through, but it is a trap for an abstract base class meant only to be subclassed, where each subclass supplies its own _run. Declare such a base with isabstract=True to opt out. See
example in the TEA tutorial.
Putting them together with simulate¶
When you call unit.simulate(), the unit runs in three stages: first _run (mass/energy balance and outlet streams), then _summary which invokes _design (size/geometry) and then _cost (purchase and installed cost). The same applies to system.simulate(), which dispatches simulate() to each unit in convergence order.
A common pitfall: calling _run directly (for debugging, or by reaching into a unit programmatically) leaves _design and _cost untouched. The unit’s outlet streams update, but its design table and cost results stay stale from the previous full simulation, which then propagates into TEA results. Always go through simulate() (on the unit or on the enclosing system) when you want results that downstream analysis will read.
Heads up. Some other tutorial cells call simulate() without explaining what it glues together; this is what they mean. See Modeling Notes & Pitfalls §4.3 for a worked symptom (a parameter sweep where cost never changes).
1.2. Useful attributes¶
Some of the class attributes that you will find useful in making your subclasses:
_N_insand_N_outsset the number of influents and effluents of theSanUnitIf you are unsure of how many influents and/or effluents there will be (e.g., they can be dynamic for a mixer), you can instead set
_ins_size_is_fixedand/or_outs_size_is_fixedtoFalse
constructionandtransportationare tuple ofConstructionandTransportationobjects for life cycle assessment (will be covered in later tutorials)purchase_cost(float) andpurchase_costs(dict) contain the total (without thes) and itemized purchase costs of this unit (i.e.,purchase_costis the sum of all the values in thepurchase_costsdict). Purchase cost of an item is calculated by multiplying the value in thebaseline_purchase_costdict with the corresponding values in theF_<D/P/M>dictSimilarly,
installed_cost(float) andinstalled_costs(dict) are the total and itemized installed costs. For each item the installed cost isbaseline_purchase_cost × (F_BM + F_D × F_P × F_M - 1), so the bare-module factorF_BMadds installation, piping, etc. on top of the design/pressure/material adjustments
add_OPEX(float ordict, USD/hr) for operating costs beyond utilities (e.g., chemicals, labor),uptime_ratio(0–1) for the fraction of time the unit runs, andlifetime(years), which drives replacement costs in TEA/LCAF_mass_in(float),mass_in(np.array containing the mass of each component),z_mass_in(np.array containing the mass fraction of each component)Additionally, there are
F_mass_out,mass_out, andz_mass_out, and the corresponding sets for molar and volume flows (e.g.,F_mol_in/F_vol_in, etc.)
H_in/H_out(changes with T) andHf_in/Hf_out(doesn’t change with T), enthalpy and enthalpy of formation of influents and effluents, respectivelyThere is also another attribute
Hnetcalculated as(H_out-H_in)+(Hf_out-Hf_in)
_graphics, how the unit will be represented when you call thediagramfunction. If not given, it will be defaulted to a boxNote that if you make the subclasses of
Mixer/Splitter/HXutility/HXprocess, the default graphics will be different because these units have their corresponding graphics
results, a method (i.e., you need to call it by<SanUnit>.results()instead of just<SanUnit>.result) to give you a quick summary of the unit design and cost results
2. Making a simple AerobicReactor¶
Alright, all those descriptions are abstract enough, and there are many details that will be best covered in an example. So let’s assume we want to design a very simple aerobic reactor that will convert 90% of the influent organics into CO2 and H2O.
[2]:
# By convention, class names use CapitalizedWords (CamelCase).
# This is the simplest possible AerobicReactor:
class AerobicReactor1(qs.SanUnit):
def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream'):
qs.SanUnit.__init__(self, ID, ins, outs, thermo, init_with)
def _run(self):
pass
def _design(self):
pass
def _cost(self):
pass
[3]:
# The default components are indeed useful!
cmps_default = qs.Components.load_default()
kwargs = {'particle_size': 'Dissolved gas',
'degradability': 'Undegradable',
'organic': False}
O2 = qs.Component('O2', search_ID='O2', **kwargs)
CO2 = qs.Component('CO2', search_ID='CO2', **kwargs)
cmps = qs.Components([*cmps_default, O2, CO2])
cmps.compile(ignore_inaccurate_molar_weight=True) # some components have inaccurate molar weight, but we don't care about that here
qs.set_thermo(cmps)
cmps.show()
CompiledComponents([
S_H2, S_CH4, S_CH3OH, S_Ac,
S_Prop, S_F, S_U_Inf, S_U_E,
C_B_Subst, C_B_BAP, C_B_UAP, C_U_Inf,
X_B_Subst, X_OHO_PHA, X_GAO_PHA, X_PAO_PHA,
X_GAO_Gly, X_PAO_Gly, X_OHO, X_AOO,
X_NOO, X_AMO, X_PAO, X_MEOLO,
X_FO, X_ACO, X_HMO, X_PRO,
X_U_Inf, X_U_OHO_E, X_U_PAO_E, X_Ig_ISS,
X_MgCO3, X_CaCO3, X_MAP, X_HAP,
X_HDP, X_FePO4, X_AlPO4, X_AlOH,
X_FeOH, X_PAO_PP_Lo, X_PAO_PP_Hi, S_NH4,
S_NO2, S_NO3, S_PO4, S_K,
S_Ca, S_Mg, S_CO3, S_N2,
S_O2, S_CAT, S_AN, H2O,
O2, CO2,
])
[4]:
# Now make a fake waste stream with these components
ws = qs.WasteStream('ws', H2O=1000, S_CH3OH=0.5, units='kg/hr')
ws.show()
WasteStream: ws
phase: 'l', T: 298.15 K, P: 101325 Pa
flow (g/hr): S_CH3OH 500
H2O 1e+06
WasteStream-specific properties:
pH : 7.0
Alkalinity : 2.5 mmol/L
COD : 498.2 mg/L
BOD : 357.2 mg/L
TC : 124.7 mg/L
TOC : 124.7 mg/L
Component concentrations (mg/L):
S_CH3OH 498.2
H2O 996432.8
[5]:
R1 = AerobicReactor1('R1', ins=ws)
R1.simulate()
R1.show()
AerobicReactor1: R1
ins...
[0] ws
phase: 'l', T: 298.15 K, P: 101325 Pa
flow (g/hr): S_CH3OH 500
H2O 1e+06
WasteStream-specific properties:
pH : 7.0
Alkalinity : 2.5 mmol/L
COD : 498.2 mg/L
BOD : 357.2 mg/L
TC : 124.7 mg/L
TOC : 124.7 mg/L
outs...
[0] ws1
phase: 'l', T: 298.15 K, P: 101325 Pa
flow: 0
WasteStream-specific properties: None for empty waste streams
OK, with these simple setups, we can “sort of” see something, but without the methods above we aren’t really doing anything useful, so let’s try to implement those methods
[6]:
class AerobicReactor2(qs.SanUnit):
def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream',
conversion=0.9, # default conversion to be 0.9,
aeration_rate=5, # assume we need 5 g/L of O2 pumped into the system
HRT=5, # hydraulic residence time being 5 hours
):
# Some standard codes you need to include for all subclasses of `SanUnit`
qs.SanUnit.__init__(self, ID, ins, outs, thermo, init_with)
# These are the unique attribures of `AerobicReactor`
self.conversion = conversion
self.aeration_rate = aeration_rate
self.HRT = HRT
# Assume a bare module factor of 2
self.F_BM = {'Tank': 2}
# Assume we'll have two influents - the waste stream and O2,
# as well as two effluents - treated waste stream and the generated CO2
_N_ins = 2
_N_outs = 2
def _run(self):
# This is equivalent to
# inf=self.ins[0]
# o2=self.ins[1]
inf, o2 = self.ins
eff, co2 = self.outs
o2.phase = co2.phase = 'g'
# Firstly let's calculate how much O2 we need,
# g/L (kg/m3) * m3/hr = kg/hr
o2_needed = self.aeration_rate * self.F_vol_in
o2.imass['O2'] = o2_needed # `imass` in kg/hr
# Mix the influent streams
eff.mix_from(self.ins)
# O2 gas turned into dissolved O2
eff.imass['S_O2'] = eff.imass['O2']
eff.imass['O2'] = 0
# Then we will want convert the organics,
# for demo purpose let's make it very simple,
# assume that we know ahead of time that
# we will only have `S_CH3OH`
# so reaction will be
# CH3OH + 1.5 O2 -> CO2 + 2H2O
# with the conversion defined by the user
x = self.conversion
converted_meoh = x * inf.imol['S_CH3OH']
consumed_o2 = 1.5 * converted_meoh
generated_co2 = converted_meoh
generated_h2o = 2 * converted_meoh
eff.imol['S_CH3OH'] -= converted_meoh
eff.imol['S_O2'] -= consumed_o2
eff.imol['H2O'] += generated_h2o
co2.imol['CO2'] = generated_co2
# Assume 5 wt% of MeOH is turned into biomass
eff.imass['X_OHO'] = 0.05 * inf.imass['S_CH3OH']
# We can (or seems more straightfoward to) move this into
# the `_design` method, but since these units won't change
# putting it here will save some simulation time
_units = {
'Volume': 'm3',
'Diameter': 'm',
'Height': 'm',
'Stainless steel': 'kg'
}
# As for the design, let's assume we will have a
# cylinder with a height-to-diameter ratio of 2:1
def _design(self):
D = self.design_results
tot_vol = self.outs[0].F_vol*self.HRT
rx_vol = tot_vol / 0.8 # assume 80% working volume
# You can certainly do `import math; math.pi`
dia = (2*rx_vol/3.14)**(1/3)
D['Volume'] = rx_vol
D['Diameter'] = dia
D['Height'] = H = 2 * dia
# Assume the tank has a thickness of 3 cm,
# we'll need the cover, but not the bottom
ss = 3.14*(dia**2)*H + 3.14/4*(dia**2)
# Assume the density is 7,500 kg/m3
D['Stainless steel'] = ss * 7500
# Let's assume that the reactor is
# made of stainless steel with a price of $3/kg
def _cost(self):
self.baseline_purchase_costs['Tank'] = \
3 * self.design_results['Stainless steel']
# Assume the electricity usage is proportional to the
# volumetric flow rate
self.power_utility.consumption = 0.1 * self.outs[0].F_vol
# Now it's a proper use of property,
# see the text enclosed in the pair of triple quotes?
# That's the documentation (e.g., the helpful prompt
# that will show up when users do
# `?AerobicReactor.conversion`)
@property
def conversion(self):
'''[float] Conversion of the organic matters in this reactor.'''
return self._conversion
@conversion.setter
def conversion(self, i):
if not 0 <= i <= 1:
# Include the specific values in the error messgae
# will often help the users (and many times you) in debugging
raise AttributeError('`conversion` must be within [0, 1], '
f'the provided value {i} is outside this range.')
self._conversion = i
[7]:
# Let's set up this unit again
R2 = AerobicReactor2('R2', ins=(ws.copy(), 'o2'), outs=('eff', 'co2'))
[8]:
# Voila!
R2.simulate()
print(R2.results())
Aerobic Reactor2 Units R2
Electricity Power kW 0.101
Cost USD/hr 0.00789
Design Volume m3 6.3
Diameter m 1.59
Height m 3.18
Stainless steel kg 2.04e+05
Purchase cost Tank USD 6.12e+05
Total purchase cost USD 6.12e+05
Utility cost USD/hr 0.00789
[9]:
R2.show()
AerobicReactor2: R2
ins...
[0] ws2
phase: 'l', T: 298.15 K, P: 101325 Pa
flow (g/hr): S_CH3OH 500
H2O 1e+06
WasteStream-specific properties:
pH : 7.0
Alkalinity : 2.5 mmol/L
COD : 498.2 mg/L
BOD : 357.2 mg/L
TC : 124.7 mg/L
TOC : 124.7 mg/L
[1] o2
phase: 'g', T: 298.15 K, P: 101325 Pa
flow: 5.02e+03 g/hr O2
WasteStream-specific properties: None for non-liquid waste streams
outs...
[0] eff
phase: 'l', T: 296.83 K, P: 101325 Pa
flow (g/hr): S_CH3OH 50
X_OHO 25
S_O2 4.34e+03
H2O 1e+06
WasteStream-specific properties:
pH : 33.5
Alkalinity : 2.5 mmol/L
COD : 74.4 mg/L
BOD : 49.6 mg/L
TC : 21.5 mg/L
TOC : 21.5 mg/L
TN : 1.7 mg/L
TP : 0.5 mg/L
TK : 0.1 mg/L
TSS : 19.3 mg/L
[1] co2
phase: 'g', T: 298.15 K, P: 101325 Pa
flow: 618 g/hr CO2
WasteStream-specific properties: None for non-liquid waste streams
[10]:
# Trying to put an unrealistic value will show our helpful message
R2.conversion = 1.1
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[10], line 2
1 # Trying to put an unrealistic value will show our helpful message
----> 2 R2.conversion = 1.1
Cell In[6], line 116, in AerobicReactor2.conversion(self, i)
112 def conversion(self, i):
113 if not 0 <= i <= 1:
114 # Include the specific values in the error messgae
115 # will often help the users (and many times you) in debugging
--> 116 raise AttributeError('`conversion` must be within [0, 1], '
117 f'the provided value {i} is outside this range.')
118 self._conversion = i
AttributeError: `conversion` must be within [0, 1], the provided value 1.1 is outside this range.
3. Other convenient features¶
We’ve done a good job in making the AerobicReactor class, but there are many helpful features that will make our lives much easier
3.1. Reactions¶
Here’s what we did for the reaction of MeOH and O2 to CO2 and H2O:
converted_meoh = x * inf.imol['S_CH3OH']
consumed_o2 = 1.5 * converted_meoh
generated_co2 = converted_meoh
generated_h2o = 2 * converted_meoh
eff.imol['S_CH3OH'] -= converted_meoh
eff.imol['S_O2'] -= consumed_o2
eff.imol['H2O'] += generated_h2o
co2.imol['CO2'] = generated_co2
For reactions like this, we can actually use Reaction to do it in a much more convenient way:
from qsdsan import Reaction as Rxn
# reaction definition reactant conversion
aerobic_rxn = Rxn('S_CH3OH + 1.5 O2 -> CO2 + 2 H2O', 'S_CH3OH', self.conversion)
If we have multiple reactions, we can use qs.ParallelReaction (if all reactions happen at once) or qs.SeriesReaction (if these reactions happen in sequence), and we can use qs.ReactionSystem to compile multiple qs.Reaction, qs.ParallelReaction, and qs.SeriesReaction together.
These classes are inherited from BioSTEAM. For more detailed instructions, refer to BioSTEAM’s documentation.
3.2. cost decorator¶
If we want to scale the cost of some equipment base on certain variables (e.g., scale the capital cost and electricity of a pump based on the flow rate, we can use the cost decorator (usage of decorator starts with the @ symbol, again recall the property decorator).
For the demo purpose, let’s assume that we need a pump for the aeration that will be scaled based on the mass flow rate of needed O2
from qsdsan.utils import cost
@cost('O2 flow rate', # the variable that the equipment is scaled on
'O2 pump', # name of the equipment
CEPCI=522, # reference-year cost index (qsdsan accepts CEPCI=; BioSTEAM calls it CE)
S=40000, # value of the scaling basis
cost=22500, # cost of the equipment when the variable is at the basis value
n=0.8, # exponential scaling factor
kW=75, # electricity usage
N=1, # number of this equipment (will be defaulted to 1 if not given)
BM=2.3) # bare module
class AerobicReactor(qs.SanUnit):
...
# Note that in the unit, you'll need to define what 'O2 flow rate' is
# in the `design_results` dict and its unit in the `_units` dict
_units = {
...,
'O2 flow rate': 'kg/hr',
}
def _design(self):
...
self.design_results['O2 flow rate'] = self.outs[1].F_mass
The scaling equations are (ub is the upper bound):
Python Aside: decorators and decorator factories (click to expand)
A decorator is a callable that takes a function (or class) and returns a (usually modified) replacement. The @deco syntax above a definition is shorthand:
@deco
def f():
...
# is equivalent to:
def f():
...
f = deco(f)
Decorators you have already met: @property (turns a method into a getter), @classmethod, @staticmethod.
Decorator factories are decorators that take arguments. @cost('O2 flow rate', ...) is not the decorator itself — it is a call that returns the decorator. The two-step desugaring is:
@cost('O2 flow rate', 'O2 pump', CEPCI=522, ...)
class AerobicReactor3(qs.SanUnit):
...
# is equivalent to:
class AerobicReactor3(qs.SanUnit):
...
the_decorator = cost('O2 flow rate', 'O2 pump', CEPCI=522, ...) # call returns the decorator
AerobicReactor3 = the_decorator(AerobicReactor3) # apply it
So when you see @something(...) with parentheses, the parentheses signal a factory; the result of that call is what wraps your function or class. You will meet several more QSDsan-specific decorator factories later — for example @process.dynamic_parameter(symbol='y_B', params={...}) and @process.kinetics(parameters={...}) in the Process tutorial, and @unit.add_specification(run=True) later in this tutorial.
Note: You can actually add the unit in the @ expression, e.g.,
@cost('O2 flow rate', ..., units='kg/hr')
But if later you define _units in the class definition by using
_units = {
'Reactor volume': 'm3',
'Diameter': 'm',
'Height': 'm',
'Stainless steel': 'kg'
}
You’ll throw away the previous definition.
What is CEPCI? The Chemical Engineering Plant Cost Index (CEPCI) tracks how the cost of building chemical process equipment changes over time (equipment, materials, labor, and so on). Cost correlations like the one above are anchored to a specific reference year, recorded by the CEPCI argument (here CEPCI=522). To restate a cost in another year’s dollars, it is scaled by \(CEPCI_{new}/CEPCI\), where \(CEPCI_{new}\) is the index for the target year (the
\(\frac{CEPCI_{new}}{CEPCI}\) term in the equation above).
The live \(CEPCI_{new}\) used for this scaling is qs.CEPCI; set it to put your cost estimates in a chosen year’s dollars. QSDsan provides the published index values by year in qs.CEPCI_by_year, so you can do, e.g., qs.CEPCI = qs.CEPCI_by_year[2023].
Note. BioSTEAM abbreviates this index as CE (the @cost argument). For consistency with qs.CEPCI/qs.CEPCI_by_year, the @cost decorator imported from qsdsan.utils also accepts CEPCI= as an alias, so @cost(..., CEPCI=522) is equivalent to @cost(..., CE=522). Both work; we use the CEPCI spelling here.
[11]:
# The current index used for cost scaling
qs.CEPCI
[11]:
567.5
[12]:
# If you want to look up CEPCI by year
qs.CEPCI_by_year
[12]:
{1980: 261,
1981: 297,
1982: 314,
1983: 317,
1984: 323,
1985: 325,
1986: 318,
1987: 324,
1988: 343,
1989: 355,
1990: 357.6,
1991: 361.3,
1992: 358.2,
1993: 359.2,
1994: 368.1,
1995: 381.1,
1996: 381.7,
1997: 386.5,
1998: 389.5,
1999: 390.6,
2000: 394.1,
2001: 394.3,
2002: 395.6,
2003: 402.0,
2004: 444.2,
2005: 468.2,
2006: 499.6,
2007: 525.4,
2008: 575.4,
2009: 521.9,
2010: 550.8,
2011: 585.7,
2012: 584.6,
2013: 567.3,
2014: 576.1,
2015: 556.8,
2016: 541.7,
2017: 567.5,
2018: 603.1,
2019: 607.5,
2020: 596.2,
2021: 708.8,
2022: 816.0,
2023: 797.9}
[ ]:
# Let's use `Reaction` and `cost` decorator here
from qsdsan import Reaction as Rxn
from qsdsan.utils import cost
@cost('O2 flow rate', 'O2 pump', CEPCI=522,
S=40000, cost=22500, n=0.8, kW=75, BM=2.3)
class AerobicReactor3(qs.SanUnit):
def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream',
conversion=0.9, aeration_rate=3, HRT=5):
qs.SanUnit.__init__(self, ID, ins, outs, thermo, init_with)
self.conversion = conversion
self.aeration_rate = aeration_rate
self.HRT = HRT
self.F_BM = {'Tank': 2}
self.aerobic_rxn = Rxn('S_CH3OH + 1.5 S_O2 -> CO2 + 2 H2O', 'S_CH3OH', conversion)
_N_ins = 2
_N_outs = 2
def _run(self):
inf, o2 = self.ins
eff, co2 = self.outs
o2.phase = co2.phase = 'g'
o2_needed = self.aeration_rate * self.F_vol_in
o2.imass['O2'] = o2_needed
eff.mix_from(self.ins)
eff.imass['S_O2'] = eff.imass['O2']
eff.imass['O2'] = 0
# Sync the reaction's conversion with the (possibly updated) `conversion`
# attribute so changing it later (e.g., from a specification) takes effect
self.aerobic_rxn.X = self.conversion
self.aerobic_rxn(eff.mol)
eff.imass['X_OHO'] = 0.05 * inf.imass['S_CH3OH']
eff.imass['S_CH3OH'] -= eff.imass['X_OHO']
_units = {
'Volume': 'm3',
'Diameter': 'm',
'Height': 'm',
'Stainless steel': 'kg',
'O2 flow rate': 'kg/hr'
}
def _design(self):
D = self.design_results
tot_vol = self.outs[0].F_vol*self.HRT
rx_vol = tot_vol / 0.8
dia = (2*rx_vol/3.14)**(1/3)
D['Volume'] = rx_vol
D['Diameter'] = dia
D['Height'] = H = 2 * dia
ss = 3.14*(dia**2)*H + 3.14/4*(dia**2)
D['Stainless steel'] = ss * 7500
D['O2 flow rate'] = self.outs[1].F_mass
def _cost(self):
self.baseline_purchase_costs['Tank'] = \
3 * self.design_results['Stainless steel']
self.power_utility.consumption = 0.1 * self.outs[0].F_vol
@property
def conversion(self):
'''[float] Conversion of the organic matters in this reactor.'''
return self._conversion
@conversion.setter
def conversion(self, i):
if not 0 <= i <= 1:
raise AttributeError('`conversion` must be within [0, 1], '
f'the provided value {i} is outside this range.')
self._conversion = i
[14]:
# Check out the results again
R3 = AerobicReactor3('R3', ins=(ws.copy(), 'o2_3'), outs=('eff_3', 'co2_3'))
R3.simulate()
print(R3.results())
Aerobic Reactor3 Units R3
Electricity Power kW 0.101
Cost USD/hr 0.00787
Design Volume m3 6.29
Diameter m 1.59
Height m 3.18
Stainless steel kg 2.04e+05
O2 flow rate kg/hr 0
Purchase cost Tank USD 6.11e+05
Total purchase cost USD 6.11e+05
Utility cost USD/hr 0.00787
[15]:
R3.show()
AerobicReactor3: R3
ins...
[0] ws3
phase: 'l', T: 298.15 K, P: 101325 Pa
flow (g/hr): S_CH3OH 500
H2O 1e+06
WasteStream-specific properties:
pH : 7.0
Alkalinity : 2.5 mmol/L
COD : 498.2 mg/L
BOD : 357.2 mg/L
TC : 124.7 mg/L
TOC : 124.7 mg/L
[1] o2_3
phase: 'g', T: 298.15 K, P: 101325 Pa
flow: 3.01e+03 g/hr O2
WasteStream-specific properties: None for non-liquid waste streams
outs...
[0] eff_3
phase: 'l', T: 297.35 K, P: 101325 Pa
flow (g/hr): S_CH3OH 25
X_OHO 25
S_O2 2.34e+03
H2O 1e+06
CO2 618
WasteStream-specific properties:
pH : 23.0
Alkalinity : 2.5 mmol/L
COD : 49.7 mg/L
BOD : 31.9 mg/L
TC : 182.8 mg/L
TOC : 15.3 mg/L
TN : 1.7 mg/L
TP : 0.5 mg/L
TK : 0.1 mg/L
TSS : 19.3 mg/L
[1] co2_3
phase: 'g', T: 298.15 K, P: 101325 Pa
flow: 0
WasteStream-specific properties: None for non-liquid waste streams
More ``cost`` arguments. A few arguments beyond the ones above are handy:
lb/ub(lower/upper size bounds) andN(number of parallel units). If the scaling basis exceedsub, the equipment is split intoN = ceil(S_new / ub)identical parallel units, each sized within bounds, rather than extrapolating one oversized unit past its valid range. (This is theNin the cost equation above.)lifetimesets a replacement interval (years) for that cost item, used by TEA/LCA.ftakes a custom cost functionf(S) -> costfor correlations that are not a simple power law; when given, it overrides thecost/npower-law form.
You can also stack multiple @cost decorators on one class: each adds its own item to baseline_purchase_costs (and its own electricity via kW), so a unit can declare a pump, a blower, and a tank as separate scaled cost items.
3.3. Auxiliary units¶
Sometimes a unit’s design includes another whole unit operation, for example a tank that needs a heat exchanger to warm its contents, or a reactor that needs a pump. Rather than re-deriving that unit’s design and cost by hand, you can make it an auxiliary unit: a full SanUnit (or BioSTEAM Unit) owned by the parent, whose design, cost, and utilities are computed by the auxiliary itself and then rolled up into the parent automatically.
To set one up: list the attribute name in the class attribute auxiliary_unit_names, build the auxiliary unit in __init__ and bind it to that name, then simulate it inside the parent’s _design. Below, a HeatedTank owns an HXutility heat exchanger that heats the influent to a target temperature.
[ ]:
from qsdsan import SanUnit
from qsdsan.unit_operations import HXutility
class HeatedTank(qs.SanUnit):
_N_ins = 1
_N_outs = 1
# declare the auxiliary slot as a class attribute
auxiliary_unit_names = ('heat_exchanger',)
def __init__(self, ID='', ins=None, outs=(), thermo=None,
init_with='WasteStream', T=320.15):
qs.SanUnit.__init__(self, ID, ins, outs, thermo, init_with)
self.T = T
self.F_BM['Tank'] = 2
# build the auxiliary heat exchanger and bind it to the declared name
# note that the name of the auxiliary unit must match the one declared in `auxiliary_unit_names`
self.heat_exchanger = HXutility(self.ID+'_hx', ins=self.ins[0].copy(), T=T)
def _run(self):
self.outs[0].copy_like(self.ins[0])
self.outs[0].T = self.T
def _design(self):
hx = self.heat_exchanger
hx.ins[0].copy_like(self.ins[0])
hx.T = self.T
hx.simulate() # the auxiliary designs and costs itself
def _cost(self):
self.baseline_purchase_costs['Tank'] = 1e4
[17]:
HT = HeatedTank('HT', ins=ws.copy())
HT.simulate()
print('auxiliary units:', [u.ID for u in HT.auxiliary_units])
print(f"tank only: {HT.baseline_purchase_costs['Tank']:,.0f} USD")
print(f"heat exchanger: {HT.heat_exchanger.purchase_cost:,.0f} USD")
print(f"unit purchase: {HT.purchase_cost:,.0f} USD (tank + auxiliary)")
print(f"unit utility: {HT.utility_cost:.3f} USD/hr (from the auxiliary's heat duty)")
auxiliary units: ['HT_hx']
tank only: 10,000 USD
heat exchanger: 3,275 USD
unit purchase: 13,275 USD (tank + auxiliary)
unit utility: 0.595 USD/hr (from the auxiliary's heat duty)
[18]:
print(HT.results())
Heated tank Units HT
Low pressure steam Duty kJ/hr 9.68e+04
Flow kmol/hr 2.5
Cost USD/hr 0.595
Purchase cost Tank USD 1e+04
Heat exchanger - Double pipe USD 3.27e+03
Total purchase cost USD 1.33e+04
Utility cost USD/hr 0.595
The auxiliary’s purchase cost is folded into the parent’s purchase_cost (and installed_cost), and its heat-exchanger steam duty shows up in the parent’s results() and utility_cost. You did not have to write any heat-exchanger design or cost code; the HXutility did it.
3.4. Equipment¶
If the same piece of equipment shows up across many of your units, factor it into an Equipment subclass. Each Equipment subclass implements its own _design and _cost; your unit then calls add_equipment_design() inside its _design and add_equipment_cost() inside its _cost to fold the equipment’s design results and costs into the unit. For a worked example, see the
ElectrochemicalCell unit and the equipment source code.
An auxiliary unit (section 3.3) and an Equipment differ in granularity: an auxiliary unit is a full owned Unit that designs and costs itself, while an Equipment is a lighter costed sub-component folded in through the add_equipment_design/add_equipment_cost hooks.
3.5. Choosing the stream class with init_with¶
Every subclass above passes init_with straight through to SanUnit.__init__. It sets which stream class each port uses ('WasteStream' by default; also 'SanStream' or 'Stream', or a dict for per-port control). See 4. SanUnit (basic) for the full string/dict forms and runnable examples.
3.6. Construction and transportation (for LCA)¶
Every SanUnit can carry construction and transportation inventories that feed life cycle assessment: Construction and Transportation objects that record quantities (e.g., kg of concrete, km of truck transport). For a custom unit that always needs certain materials, declare them as a class attribute via _construction_specs (resolved lazily when an LCA is created), or set unit.construction / unit.transportation on the instance. The full workflow, including the
ImpactItems these link to, is covered in 8. LCA.
3.7. Specifications¶
A specification is a function attached to a unit that runs during simulation. Use it to adjust a unit, or drive it toward a target, without writing a whole subclass. Add one with add_specification (often as a decorator); pass run=True to also call the unit’s _run afterward.
[19]:
# Reuse the AerobicReactor3 class from above, but tune its conversion on the fly
U_spec = AerobicReactor3('U_spec', ins=(ws.copy(), 'o2_spec'), outs=('eff_spec', 'co2_spec'))
U_spec.simulate() # default conversion (0.9)
print(f'No spec: conversion {U_spec.conversion}, effluent COD {U_spec.outs[0].COD:.1f} mg/L')
@U_spec.add_specification(run=True) # run=True: also run the unit after the spec
def boost_conversion():
# use a higher conversion when the influent is more concentrated
U_spec.conversion = 0.95 if U_spec.ins[0].COD > 400 else 0.85
U_spec.simulate()
print(f'With spec: conversion {U_spec.conversion}, effluent COD {U_spec.outs[0].COD:.1f} mg/L')
No spec: conversion 0.9, effluent COD 49.7 mg/L
With spec: conversion 0.95, effluent COD 24.6 mg/L
The spec sets a higher conversion for the concentrated influent, and the effluent COD drops accordingly between the two runs above. The run=True flag makes sure the change is implemented for this round of simulation: after the spec function sets conversion, the unit’s _run is called again, so the new value propagates to the outlets (and then to _design/_cost) within the same simulate(). With the default run=False the spec is still called, but the unit’s mass
balance is not re-evaluated, so a parameter the spec changed would not reach the outlets until the next _run.
Degrees of freedom. A specification is one free hand on the unit. To pin an outlet (or any result) to a target value, you free exactly one adjustable parameter for the spec to set: one parameter per target. Two independent targets need two parameters, and asking for more targets than the unit has free parameters is over-constrained, with no consistent solution. The same idea scales up to whole systems, covered in System tutorial section 3.
↑ Back to top