diff --git a/exposan/metab/README.rst b/exposan/metab/README.rst new file mode 100644 index 00000000..43f9c3c2 --- /dev/null +++ b/exposan/metab/README.rst @@ -0,0 +1,117 @@ +================================================================= +metab: Modular Encapsulated Two-stage Anaerobic Biological system +================================================================= + +Summary +------- +This module serves as a process simulator for the Modular Encapsulated Two-stage Anaerobic Biological (METAB) system and its design variants, +as described in Zhang et al. [1]_ + +This module include a "flexible" version of the Anaerobic Digestion Model no.1 (ADM1) [2]_, which enables ideal control of pH and estimation +of the required acid/base dosing by directly setting a target pH. + +Multiple reactor types for the METAB system are included in the ``units.py`` file in this module: + + - upflow sludge blanket reactor (``UASB``) + - fluidized bed reactor (``METAB_FluidizedBed``) + - packed bed reactor (``METAB_PackedBed``) + - batch experiment reactor (``METAB_BatchExp``) + +Among them, ``METAB_FluidizedBed``, ``METAB_PackedBed``, and ``METAB_BatchExp`` include a micro-scale mechanistic model of mass transfer through +encapsulation matrix in reactor mass balance. + +.. figure:: ./readme_figures/example_system.png + + *An example METAB system layout: two packed bed reactors in series with sidestream membrane gas extraction at the 1st stage and effluent degassing.* + +To reproduce the results and figures included in Zhang et al. [1]_, run corresponding functions in ``analyses.py``. You can find a full list of the +packages in the environment used to generate the results in `qsdsan.yml `_. + +Getting Started +--------------- +.. code-block:: python + + >>> # create a METAB system by specifying reactor type, gas extraction method, design HRT etc. + >>> from exposan.metab import create_system + >>> sys = create_system( + ... n_stages=2, # number of stages + ... reactor_type='PB', # PB for packed bed, FB for fluidized bed, or UASB + ... gas_extraction='M', # M for membrane gas extraction, V for vacuum extraction, P for passive venting + ... Q=5, # influent flowrate in m3/d + ... T=22, # reactor temperature in degree C + ... tot_HRT=12, # total HRT in d + ... ) + + >>> # Simulate the system for 200 days + >>> sys.simulate(state_reset_hook='reset_cache', t_span=(0,200), method='BDF') + System: PB2M_edg + Highest convergence error among components in recycle + streams {DMs-1} after 3 loops: + - flow rate 1.82e-12 kmol/hr (4e-14%) + - temperature 0.00e+00 K (0%) + ins... + [0] inf + phase: 'l', T: 295.15 K, P: 101325 Pa + flow (kmol/hr): S_su 0.00347 + S_aa 0.125 + S_fa 0.000325 + S_va 0.000816 + S_bu 0.000957 + S_pro 0.00114 + S_ac 0.00141 + ... 11.6 + outs... + [0] bg1 + phase: 'g', T: 295.15 K, P: 101325 Pa + flow: 0 + [1] bg2 + phase: 'g', T: 295.15 K, P: 101325 Pa + flow (kmol/hr): S_h2 1.49e-07 + S_ch4 0.0589 + S_IC 0.000946 + H2O 0.000489 + [2] bge + phase: 'g', T: 295.15 K, P: 101325 Pa + flow (kmol/hr): S_h2 3.76e-10 + S_ch4 0.000449 + S_IC 0.000413 + [3] eff_dg + phase: 'l', T: 295.15 K, P: 101325 Pa + flow (kmol/hr): S_su 2.28e-07 + S_aa 4.9e-05 + S_fa 7.64e-06 + S_va 5.64e-07 + S_bu 8.55e-07 + S_pro 8.7e-07 + S_ac 7.4e-06 + ... 11.7 + [4] bgs + phase: 'g', T: 295.15 K, P: 101325 Pa + flow (kmol/hr): S_h2 1.01e-05 + S_ch4 0.0213 + S_IC 0.00291 + + >>> # You can look at TEA and LCA results after simulation + >>> sys.TEA + TEA: PB2M_edg + NPV : -1,413,755 USD at 10.0% discount rate + + >>> sys.LCA + LCA: PB2M_edg (lifetime 30 yr) + Impacts: + Construction Transportation Stream Others Total + ODP (kg CFC-11-Eq) 0.0196 0 -0.0322 0.0252 0.0126 + HTNC (CTUh) 0.18 0 -0.0046 0.121 0.296 + EP (kg N-Eq) 1.65e+03 0 -27.6 3.24e+03 4.87e+03 + PMFP (PM2.5-Eq) 467 0 -16.3 1.27e+03 1.72e+03 + MIR (kg O3-Eq) 2.92e+04 0 -2.32e+03 1.4e+04 4.09e+04 + EF (CTUe) 1.25e+07 0 -2.33e+05 5.11e+06 1.74e+07 + GWP100 (kg CO2-Eq) 5.6e+05 0 -3.78e+04 4.45e+05 9.67e+05 + HTC (CTUh) 0.061 0 -0.00259 0.0375 0.0959 + AP (kg SO2-Eq) 2.52e+03 0 -205 946 3.26e+03 + + +References +---------- +.. [1] Zhang et al., Pathway toward sustainable distributed treatment of high strength food industry wastewater with encapsulated anaerobic technology. 2023, *In preparation* +.. [2] IWA Task Group for Mathematical Modelling of Anaerobic Digestion Processes. Anaerobic Digestion Model No.1 (ADM1); IWA Publishing, 2005. ``_ diff --git a/exposan/metab/qsdsan.yml b/exposan/metab/qsdsan.yml new file mode 100644 index 00000000..bc3f9bce --- /dev/null +++ b/exposan/metab/qsdsan.yml @@ -0,0 +1,296 @@ +name: dev +channels: + - defaults +dependencies: + - alabaster=0.7.12=pyhd3eb1b0_0 + - arrow=1.2.3=py39haa95532_1 + - astroid=2.14.2=py39haa95532_0 + - asttokens=2.0.5=pyhd3eb1b0_0 + - atomicwrites=1.4.0=py_0 + - attrs=22.1.0=py39haa95532_0 + - autopep8=1.6.0=pyhd3eb1b0_1 + - babel=2.11.0=py39haa95532_0 + - backcall=0.2.0=pyhd3eb1b0_0 + - bcrypt=3.2.0=py39h2bbff1b_1 + - beautifulsoup4=4.12.2=py39haa95532_0 + - binaryornot=0.4.4=pyhd3eb1b0_1 + - black=23.3.0=py39haa95532_0 + - bleach=4.1.0=pyhd3eb1b0_0 + - boost-cpp=1.73.0=h2bbff1b_12 + - brotlipy=0.7.0=py39h2bbff1b_1003 + - ca-certificates=2023.01.10=haa95532_0 + - cairo=1.16.0=haedb8bc_4 + - certifi=2022.12.7=py39haa95532_0 + - cffi=1.15.1=py39h2bbff1b_3 + - chardet=4.0.0=py39haa95532_1003 + - charset-normalizer=2.0.4=pyhd3eb1b0_0 + - click=8.0.4=py39haa95532_0 + - cloudpickle=2.0.0=pyhd3eb1b0_0 + - colorama=0.4.6=py39haa95532_0 + - cookiecutter=1.7.3=pyhd3eb1b0_0 + - cryptography=39.0.1=py39h21b164f_0 + - debugpy=1.5.1=py39hd77b12b_0 + - decorator=5.1.1=pyhd3eb1b0_0 + - defusedxml=0.7.1=pyhd3eb1b0_0 + - diff-match-patch=20200713=pyhd3eb1b0_0 + - dill=0.3.6=py39haa95532_0 + - docstring-to-markdown=0.11=py39haa95532_0 + - docutils=0.18.1=py39haa95532_3 + - entrypoints=0.4=py39haa95532_0 + - executing=0.8.3=pyhd3eb1b0_0 + - expat=2.4.9=h6c2663c_0 + - flake8=6.0.0=py39haa95532_0 + - font-ttf-dejavu-sans-mono=2.37=hd3eb1b0_0 + - font-ttf-inconsolata=2.001=hcb22688_0 + - font-ttf-source-code-pro=2.030=hd3eb1b0_0 + - font-ttf-ubuntu=0.83=h8b1ccd4_0 + - fontconfig=2.14.1=h9c4af85_2 + - fonts-anaconda=1=h8fa9717_0 + - fonts-conda-ecosystem=1=hd3eb1b0_0 + - freetype=2.12.1=ha860e81_0 + - fribidi=1.0.10=h62dcd97_0 + - getopt-win32=0.1=h2bbff1b_0 + - giflib=5.2.1=h8cc25b3_3 + - glib=2.69.1=h5dc1a3c_2 + - graphite2=1.3.14=hd77b12b_1 + - graphviz=2.50.0=h7eca76f_1 + - gst-plugins-base=1.18.5=h9e645db_0 + - gstreamer=1.18.5=hd78058f_0 + - gts=0.7.6=h63ab5a1_3 + - harfbuzz=4.3.0=hb646838_1 + - icu=58.2=ha925a31_3 + - idna=3.4=py39haa95532_0 + - imagesize=1.4.1=py39haa95532_0 + - importlib-metadata=6.0.0=py39haa95532_0 + - importlib_metadata=6.0.0=hd3eb1b0_0 + - inflection=0.5.1=py39haa95532_0 + - intervaltree=3.1.0=pyhd3eb1b0_0 + - ipykernel=6.19.2=py39hd4e2768_0 + - ipython=8.12.0=py39haa95532_0 + - ipython_genutils=0.2.0=pyhd3eb1b0_1 + - isort=5.9.3=pyhd3eb1b0_0 + - jaraco.classes=3.2.1=pyhd3eb1b0_0 + - jedi=0.18.1=py39haa95532_1 + - jellyfish=0.9.0=py39h2bbff1b_0 + - jinja2=3.1.2=py39haa95532_0 + - jinja2-time=0.2.0=pyhd3eb1b0_3 + - jpeg=9e=h2bbff1b_1 + - jsonschema=4.17.3=py39haa95532_0 + - jupyter_client=8.1.0=py39haa95532_0 + - jupyter_core=5.3.0=py39haa95532_0 + - jupyterlab_pygments=0.1.2=py_0 + - keyring=23.13.1=py39haa95532_0 + - krb5=1.19.4=h5b6d351_0 + - lazy-object-proxy=1.6.0=py39h2bbff1b_0 + - lcms2=2.12=h83e58a3_0 + - lerc=3.0=hd77b12b_0 + - libboost=1.73.0=h6c2663c_12 + - libclang=14.0.6=default_hb5a9fac_1 + - libclang13=14.0.6=default_h8e68704_1 + - libcurl=7.88.1=h86230a5_0 + - libdeflate=1.17=h2bbff1b_0 + - libffi=3.4.2=hd77b12b_6 + - libgd=2.3.3=hd77b12b_2 + - libiconv=1.16=h2bbff1b_2 + - libogg=1.3.5=h2bbff1b_1 + - libpng=1.6.39=h8cc25b3_0 + - libsodium=1.0.18=h62dcd97_0 + - libspatialindex=1.9.3=h6c2663c_0 + - libssh2=1.10.0=hcd4344a_0 + - libtiff=4.5.0=h6c2663c_2 + - libvorbis=1.3.7=he774522_0 + - libwebp=1.2.4=hbc33d0d_1 + - libwebp-base=1.2.4=h2bbff1b_1 + - libxml2=2.10.3=h0ad7f3c_0 + - libxslt=1.1.37=h2bbff1b_0 + - lxml=4.9.2=py39h2bbff1b_0 + - lz4-c=1.9.4=h2bbff1b_0 + - markupsafe=2.1.1=py39h2bbff1b_0 + - matplotlib-inline=0.1.6=py39haa95532_0 + - mccabe=0.7.0=pyhd3eb1b0_0 + - mistune=0.8.4=py39h2bbff1b_1000 + - more-itertools=8.12.0=pyhd3eb1b0_0 + - mypy_extensions=0.4.3=py39haa95532_1 + - nbclient=0.5.13=py39haa95532_0 + - nbconvert=6.5.4=py39haa95532_0 + - nbformat=5.7.0=py39haa95532_0 + - nest-asyncio=1.5.6=py39haa95532_0 + - numpydoc=1.5.0=py39haa95532_0 + - openjpeg=2.4.0=h4fc8c34_0 + - openssl=1.1.1t=h2bbff1b_0 + - packaging=23.0=py39haa95532_0 + - pandocfilters=1.5.0=pyhd3eb1b0_0 + - pango=1.50.7=h78c2152_0 + - paramiko=2.8.1=pyhd3eb1b0_0 + - parso=0.8.3=pyhd3eb1b0_0 + - pathspec=0.10.3=py39haa95532_0 + - pcre=8.45=hd77b12b_0 + - pexpect=4.8.0=pyhd3eb1b0_3 + - pickleshare=0.7.5=pyhd3eb1b0_1003 + - pip=23.0.1=py39haa95532_0 + - pixman=0.40.0=h2bbff1b_1 + - platformdirs=2.5.2=py39haa95532_0 + - pluggy=1.0.0=py39haa95532_1 + - ply=3.11=py39haa95532_0 + - poppler=22.12.0=h268424c_0 + - poppler-data=0.4.11=haa95532_1 + - poyo=0.5.0=pyhd3eb1b0_0 + - prompt-toolkit=3.0.36=py39haa95532_0 + - psutil=5.9.0=py39h2bbff1b_0 + - ptyprocess=0.7.0=pyhd3eb1b0_2 + - pure_eval=0.2.2=pyhd3eb1b0_0 + - pycodestyle=2.10.0=py39haa95532_0 + - pycparser=2.21=pyhd3eb1b0_0 + - pydocstyle=6.3.0=py39haa95532_0 + - pyflakes=3.0.1=py39haa95532_0 + - pylint=2.16.2=py39haa95532_0 + - pylint-venv=2.3.0=py39haa95532_0 + - pyls-spyder=0.4.0=pyhd3eb1b0_0 + - pynacl=1.5.0=py39h8cc25b3_0 + - pyopenssl=23.0.0=py39haa95532_0 + - pyqt=5.15.7=py39hd77b12b_0 + - pyqt5-sip=12.11.0=py39hd77b12b_0 + - pyqtwebengine=5.15.7=py39hd77b12b_0 + - pyrsistent=0.18.0=py39h196d8e1_0 + - pysocks=1.7.1=py39haa95532_0 + - python=3.9.16=h6244533_2 + - python-dateutil=2.8.2=pyhd3eb1b0_0 + - python-fastjsonschema=2.16.2=py39haa95532_0 + - python-lsp-black=1.2.1=py39haa95532_0 + - python-lsp-jsonrpc=1.0.0=pyhd3eb1b0_0 + - python-lsp-server=1.7.2=py39haa95532_0 + - python-slugify=5.0.2=pyhd3eb1b0_0 + - pytoolconfig=1.2.5=py39haa95532_1 + - pytz=2022.7=py39haa95532_0 + - pywin32=305=py39h2bbff1b_0 + - pywin32-ctypes=0.2.0=py39haa95532_1000 + - pyyaml=6.0=py39h2bbff1b_1 + - pyzmq=25.0.2=py39hd77b12b_0 + - qdarkstyle=3.0.2=pyhd3eb1b0_0 + - qstylizer=0.2.2=py39haa95532_0 + - qt-main=5.15.2=he8e5bd7_8 + - qt-webengine=5.15.9=hb9a9bb5_5 + - qtawesome=1.2.2=py39haa95532_0 + - qtconsole=5.4.2=py39haa95532_0 + - qtpy=2.2.0=py39haa95532_0 + - qtwebkit=5.212=h2bbfb41_5 + - requests=2.29.0=py39haa95532_0 + - rope=1.7.0=py39haa95532_0 + - rtree=1.0.1=py39h2eaa2aa_0 + - setuptools=66.0.0=py39haa95532_0 + - sip=6.6.2=py39hd77b12b_0 + - six=1.16.0=pyhd3eb1b0_1 + - snowballstemmer=2.2.0=pyhd3eb1b0_0 + - sortedcontainers=2.4.0=pyhd3eb1b0_0 + - soupsieve=2.4=py39haa95532_0 + - sphinxcontrib-applehelp=1.0.2=pyhd3eb1b0_0 + - sphinxcontrib-devhelp=1.0.2=pyhd3eb1b0_0 + - sphinxcontrib-htmlhelp=2.0.0=pyhd3eb1b0_0 + - sphinxcontrib-jsmath=1.0.1=pyhd3eb1b0_0 + - sphinxcontrib-qthelp=1.0.3=pyhd3eb1b0_0 + - spyder=5.4.3=py39haa95532_1 + - spyder-kernels=2.4.3=py39haa95532_0 + - sqlite=3.41.2=h2bbff1b_0 + - stack_data=0.2.0=pyhd3eb1b0_0 + - text-unidecode=1.3=pyhd3eb1b0_0 + - textdistance=4.2.1=pyhd3eb1b0_0 + - three-merge=0.1.1=pyhd3eb1b0_0 + - tinycss2=1.2.1=py39haa95532_0 + - toml=0.10.2=pyhd3eb1b0_0 + - tomli=2.0.1=py39haa95532_0 + - tomlkit=0.11.1=py39haa95532_0 + - tornado=6.2=py39h2bbff1b_0 + - traitlets=5.7.1=py39haa95532_0 + - typing-extensions=4.5.0=py39haa95532_0 + - typing_extensions=4.5.0=py39haa95532_0 + - ujson=5.4.0=py39hd77b12b_0 + - unidecode=1.2.0=pyhd3eb1b0_0 + - urllib3=1.26.15=py39haa95532_0 + - vc=14.2=h21ff451_1 + - vs2015_runtime=14.27.29016=h5e58377_2 + - watchdog=2.1.6=py39haa95532_0 + - wcwidth=0.2.5=pyhd3eb1b0_0 + - webencodings=0.5.1=py39haa95532_1 + - whatthepatch=1.0.2=py39haa95532_0 + - wheel=0.38.4=py39haa95532_0 + - win_inet_pton=1.1.0=py39haa95532_0 + - wrapt=1.14.1=py39h2bbff1b_0 + - xz=5.2.10=h8cc25b3_1 + - yaml=0.2.5=he774522_0 + - yapf=0.31.0=pyhd3eb1b0_0 + - zeromq=4.3.4=hd77b12b_0 + - zipp=3.11.0=py39haa95532_0 + - zlib=1.2.13=h8cc25b3_0 + - zstd=1.5.5=hd43e919_0 + - pip: + - ansicolors==1.1.8 + - chaospy==4.3.12 + - chemicals==1.1.4 + - colorpalette==0.3.3 + - comm==0.1.4 + - contourpy==1.0.7 + - coverage==7.3.1 + - cycler==0.11.0 + - cython==0.29.34 + - dcor==0.6 + - et-xmlfile==1.1.0 + - exceptiongroup==1.1.3 + - fdasrsf==2.4.0 + - findiff==0.9.2 + - flexsolve==0.5.4 + - fluids==1.0.23 + - fonttools==4.39.3 + - free-properties==0.3.6 + - furo==2023.9.10 + - gpy==1.12.0 + - importlib-resources==5.12.0 + - iniconfig==2.0.0 + - ipywidgets==8.1.1 + - joblib==1.2.0 + - jupyterlab-widgets==3.0.9 + - kiwisolver==1.4.4 + - lazy-loader==0.2 + - llvmlite==0.39.1 + - matplotlib==3.8.0 + - mpmath==1.3.0 + - multimethod==1.9.1 + - multiprocess==0.70.14 + - nbsphinx==0.9.3 + - nbval==0.10.0 + - numba==0.56.4 + - numpoly==1.2.7 + - numpy==1.23.5 + - openpyxl==3.1.2 + - pandas==2.0.1 + - pandoc==2.3 + - paramz==0.9.5 + - patsy==0.5.3 + - pillow==9.5.0 + - pint==0.20.1 + - plumbum==1.8.2 + - pyglet==2.0.5 + - pygments==2.16.1 + - pyparsing==3.0.9 + - pytest==7.4.2 + - python-graphviz==0.20.1 + - rdata==0.9 + - salib==1.4.7 + - scikit-datasets==0.2.2 + - scikit-fda==0.8.1 + - scikit-learn==1.2.2 + - scipy==1.10.1 + - seaborn==0.13.0 + - sphinx==7.2.6 + - sphinx-basic-ng==1.0.0b2 + - sphinx-copybutton==0.5.2 + - sphinx-design==0.5.0 + - sphinxcontrib-serializinghtml==1.1.9 + - sympy==1.11.1 + - thermo==0.2.26 + - threadpoolctl==3.1.0 + - tqdm==4.65.0 + - tzdata==2023.3 + - widgetsnbextension==4.0.9 + - xarray==2023.4.2 + - xlrd==2.0.1 +prefix: C:\Users\joy_c\anaconda3\envs\dev diff --git a/exposan/metab/readme_figures/example_system.png b/exposan/metab/readme_figures/example_system.png new file mode 100644 index 00000000..ba6bde54 Binary files /dev/null and b/exposan/metab/readme_figures/example_system.png differ diff --git a/exposan/metab/units.py b/exposan/metab/units.py index 8676a753..ca0b4d4b 100644 --- a/exposan/metab/units.py +++ b/exposan/metab/units.py @@ -34,7 +34,8 @@ __all__ = ('DegassingMembrane', 'UASB', 'METAB_FluidizedBed', - 'METAB_PackedBed') + 'METAB_PackedBed', + 'METAB_BatchExp') _fts2mhr = auom('ft/s').conversion_factor('m/hr') _cmph_2_gpm = auom('m3/hr').conversion_factor('gpm') @@ -455,12 +456,16 @@ def _compile_ODE(self): else: f_qgas = self.f_q_gas_var_P_headspace def dy_dt(t, QC_ins, QC, dQC_ins): + #!!! to avoid accumulation of floating error due to limited precision + QC[np.abs(QC) < 2.22044604925e-16] = 0 S_liq = QC[:n_cmps] S_gas = QC[n_cmps: (n_cmps+n_gas)] Q_ins = QC_ins[:, -1] S_ins = QC_ins[:, :-1] * 1e-3 # mg/L to kg/m3 Q = sum(Q_ins) rhos = _f_rhos(QC) + #!!! to avoid accumulation of floating error due to limited precision + rhos[np.abs(rhos) < 2.22044604925e-16] = 0 _dstate[:n_cmps] = (Q_ins @ S_ins - Q*S_liq*(1-f_rtn))/V_liq \ + np.dot(M_stoichio.T, rhos) q_gas = f_qgas(rhos[-3:], S_gas, T) @@ -1616,4 +1621,209 @@ def biomass_tss(self, biomass_IDs): C_en_avg = np.dot(en_bm, dV)/V_bead outs.append([bk_bm, C_en_avg]) outs = np.asarray(outs) - return np.mean(outs, axis=0) \ No newline at end of file + return np.mean(outs, axis=0) + +#%% Batch experiment + +class METAB_BatchExp(METAB_FluidizedBed): + _N_ins = 0 + _N_outs = 1 + + def __init__(self, ID='', outs=(), thermo=None, + init_with='WasteStream', V_liq=25, V_gas=30, + V_beads=10, bead_diameter=10, n_layer=5, + boundary_layer_thickness=0.01, diffusivity=None, + f_diff=0.55, max_encapsulation_tss=16, model=None, + pH_ctrl=False, T=295.15, headspace_P=1.013, external_P=1.013, + pipe_resistance=5.0e-1, fixed_headspace_P=True, + isdynamic=True, exogenous_vars=(), **kwargs): + + SanUnit.__init__(self, ID=ID, ins=None, outs=outs, thermo=thermo, + init_with=init_with, isdynamic=isdynamic, + exogenous_vars=exogenous_vars, **kwargs) + self.V_gas = V_gas + self.V_liq = V_liq + self.V_beads = V_beads + self.bead_diameter = bead_diameter + self.n_layer = n_layer + self.boundary_layer_thickness = boundary_layer_thickness + self.diffusivity = diffusivity + self.f_diff = f_diff + self.max_encapsulation_tss = max_encapsulation_tss + + self.pH_ctrl = pH_ctrl + self.T = T + self._q_gas = 0 + self._n_gas = None + self._gas_cmp_idx = None + self._state_keys = None + self._S_vapor = None + self._biogas = WasteStream(phase='g') + self.headspace_P = headspace_P + self.external_P = external_P + self.pipe_resistance = pipe_resistance + self.fixed_headspace_P = fixed_headspace_P + self._tempstate = [] + + self.model = model + self._cached_state = None + + @property + def V_beads(self): + return self._V_beads + @V_beads.setter + def V_beads(self, V): + self._V_beads = V + + def _setup(self): + AnaerobicCSTR._setup(self) + + def _run(self): + pass + + def _set_init_Cs(self, arr=None, **kwargs): + cmps = self.thermo.chemicals + cmpx = cmps.index + if arr is None: + Cs = np.zeros(len(cmps)) + for k, v in kwargs.items(): Cs[cmpx(k)] = v + return Cs + else: + arr = np.asarray(arr) + if arr.shape != (len(cmps),): + raise ValueError(f'arr must be None or a 1d array of length {len(cmps)}, ' + f'not {arr.shape}') + return arr + + def set_bulk_init(self, arr=None, **kwargs): + self._bulk_concs = self._set_init_Cs(arr, **kwargs) + + def set_encap_init(self, arr=None, **kwargs): + n_dz = self.n_dz + Cs = self._set_init_Cs(arr, **kwargs) + self._encap_concs = np.tile(Cs, n_dz) + + def _init_state(self): + if self._cached_state is not None: + Cs = self._cached_state + Cs[-len(self._bulk_concs):] = self._bulk_concs + else: Cs = np.append(self._encap_concs, self._bulk_concs) + self._state = np.append(Cs, [0]*self._n_gas + [0]).astype('float64') + self._dstate = self._state * 0. + + def _update_state(self): + cmps = self.components + n_cmps = len(cmps) + n_gas = self._n_gas + y_gas = self._state[-(n_gas+1):-1] + i_mass = cmps.i_mass + chem_MW = cmps.chem_MW + gas, = self._outs + if gas.state is None: + gas.state = np.zeros(n_cmps+1) + gas.state[self._gas_cmp_idx] = y_gas + gas.state[cmps.index('H2O')] = self._S_vapor + gas.state[-1] = self._q_gas + gas.state[:n_cmps] = gas.state[:n_cmps] * chem_MW / i_mass * 1e3 # i.e., M biogas to mg (measured_unit) / L + + def _update_dstate(self): + self._tempstate = self.model.rate_function._params['root'].data.copy() + gas, = self._outs + if gas.dstate is None: + # contains no info on dstate + n_cmps = len(self.components) + gas.dstate = np.zeros(n_cmps+1) + + def _compile_ODE(self): + cmps = self.components + _dstate = self._dstate + _update_dstate = self._update_dstate + n_cmps = len(cmps) + n_dz = self.n_layer + n_gas = self._n_gas + T = self.T + + _f_rhos = self.model.flex_rate_function + _params = self.model.rate_function.params + stoi_bk = self.model.stoichio_eval() + stoi_en = stoi_bk[:-n_gas] # no liquid-gas transfer + Rho_bk = lambda state_arr: _f_rhos(state_arr, _params, + T_op=T, pH=self.pH_ctrl, + gas_transfer=True) + Rho_en = lambda state_arr: _f_rhos(state_arr, _params, + T_op=T, pH=self.pH_ctrl, + gas_transfer=False) + gas_mass2mol_conversion = (cmps.i_mass / cmps.chem_MW)[self._gas_cmp_idx] + + V_liq = self.V_liq + V_gas = self.V_gas + V_beads = self.V_beads + r_beads = self.r_beads + A_beads = 3 * V_beads / r_beads # m2, total bead surface area + + dz = r_beads / n_dz + zs = np.linspace(dz, r_beads, n_dz) + dV = 4/3*np.pi*(zs)**3 + dV[1:] -= dV[:-1] + V_bead = (4/3*np.pi*r_beads**3) + + D = self.D # Diffusivity in beads + k = self.k_bl # mass transfer coeffient through liquid boundary layer + K_tss = self.K_tss + S_idx = list(*(D.nonzero())) + D_ov_dz2 = D[S_idx]/(dz**2) # (n_soluble,) + D_ov_dz = D[S_idx]/dz + _1_ov_z = 1/zs # (n_dz,) + + if self._fixed_P_gas: + f_qgas = self.f_q_gas_fixed_P_headspace + else: + f_qgas = self.f_q_gas_var_P_headspace + + def dy_dt(t, y_ins, y, dy_ins): + S_gas = y[-(n_gas+1):-1] + Cs_bk = y[-(n_cmps+n_gas+1):-(n_gas+1)] # bulk liquid concentrations + Cs_en = y[:n_dz*n_cmps].reshape((n_dz, n_cmps)) # each row is one control volume + + # Transformation + rhos_en = np.apply_along_axis(Rho_en, 1, Cs_en) + Rs_en = rhos_en @ stoi_en # n_dz * n_cmps + rhos_bk = Rho_bk(y[-(n_cmps+n_gas+1):-1]) + Rs_bk = np.dot(stoi_bk.T, rhos_bk) # n_cmps+5 + q_gas = f_qgas(rhos_bk[-n_gas:], S_gas, T) + gas_transfer = - q_gas*S_gas/V_gas + rhos_bk[-n_gas:] * V_liq/V_gas * gas_mass2mol_conversion + + # Detachment -- particulates + tss = np.sum(Cs_en * (cmps.x*cmps.i_mass), axis=1) + x_net_growth = np.sum(Rs_en * cmps.x, axis=1)/np.sum(Cs_en * cmps.x, axis=1) # d^(-1), equivalent to k_de + u_de = 1/(1+np.exp(K_tss-tss)) * np.maximum(x_net_growth, 0) * (tss > 0) + de_en = np.diag(u_de) @ (Cs_en * cmps.x) + tot_de = np.sum(np.diag(dV) @ de_en, axis=0) / V_bead # detachment per unit volume of beads + + #!!! Mass transfer (centered differences) -- MOL; solubles only + C_lf = Cs_en[-1] + J_lf = k*(Cs_bk - C_lf) + S_en = Cs_en[:, S_idx] + M_transfer = np.zeros_like(Cs_en) + M_transfer[1:-1, S_idx] = D_ov_dz2 * (S_en[2:] - 2*S_en[1:-1] + S_en[:-2])\ + + D_ov_dz * (np.diag(_1_ov_z[1:-1]) @ (S_en[2:] - S_en[:-2])) + M_transfer[0, S_idx] = 2 * D_ov_dz2 * (S_en[1] - S_en[0]) + M_transfer[-1, S_idx] = 2 * D_ov_dz2 * (S_en[-2] - S_en[-1])\ + + 2 * (1/dz + _1_ov_z[-1]) * J_lf[S_idx] + + # Mass balance + dCdt_en = M_transfer + Rs_en - de_en + dCdt_bk = - A_beads/V_liq*J_lf + Rs_bk + V_beads*tot_de/V_liq + + _dstate[:n_dz*n_cmps] = dCdt_en.flatten() + _dstate[-(n_cmps+n_gas+1):-(n_gas+1)] = dCdt_bk + _dstate[-(n_gas+1):-1] = gas_transfer + _update_dstate() + + self._ODE = dy_dt + + def _design(self): + pass + + def _cost(self): + pass \ No newline at end of file diff --git a/exposan/metab/validation.py b/exposan/metab/validation.py new file mode 100644 index 00000000..d950279c --- /dev/null +++ b/exposan/metab/validation.py @@ -0,0 +1,235 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Aug 29 09:19:24 2023 + +@author: Joy Zhang +""" + +import qsdsan as qs, matplotlib.pyplot as plt, numpy as np +from qsdsan.utils import ospath +from qsdsan import processes as pc +from exposan.metab import flex_rhos_adm1, METAB_BatchExp, figures_path + +#%% +cmps = pc.create_adm1_cmps() +adm1 = pc.ADM1(flex_rate_function=flex_rhos_adm1) +# adm1.rate_function.params['rate_constants'][1:4] /= 16 +# adm1.rate_function.params['rate_constants'][4:10] *= 1.5 +# adm1.rate_function.params['half_sat_coeffs'][:-2] *= 0.5 +# adm1.rate_function.params['rate_constants'][10:12] *= 1.5 +# adm1.rate_function.params['half_sat_coeffs'][-2:] *= 0.5 +# adm1.rate_function.params['rate_constants'][12:19] *= 2 + +# adm1.set_parameters( +# # f_bu_aa=0.18, f_pro_aa=0.35, f_ac_aa=0.33, f_va_aa=0.14 +# # 'f_va_aa': 0.23, +# # 'f_bu_aa': 0.26, +# # 'f_pro_aa': 0.05, +# # 'f_ac_aa': 0.4, +# # 'f_h2_aa': 0.06, +# # Y_aa=0.04, +# # Y_ac=0.025 +# ) + +syn_ww = dict( + # X_ch=2.065, + # X_pr=5.832, + X_ch=2.033, + X_pr=5.812, + X_li=0.004, + S_su=1.4e-3, + S_aa=0.159, + S_fa=2.755, + S_IC=0.0143, + # S_IN=0.014, + S_cat=0.0075, + S_an=0.0058 + ) + +#%% +from exposan.metab import create_system +sys = create_system(1, inf_concs=syn_ww, tot_HRT=8) +# mdl = sys.units[0].model +# mdl.rate_function.params['rate_constants'][1:4] /= 16 +# mdl.rate_function.params['rate_constants'][4:12] *= 1.5 +# mdl.rate_function.params['half_sat_coeffs'][:-2] *= 0.5 +sys.simulate(state_reset_hook='reset_cache', t_span=(0,200), method='BDF') + +#%% +# yields ~ 0.1g/L VSS +# h2r_seed = dict( +# X_su=4.69/59, +# X_aa=1.23/59, +# X_fa=0.484/59, +# X_c4=1.42/59, +# X_pro=0.898/59 +# ) + +h2r_seed = { + 'X_su': 1.956e-2/0.7, + 'X_aa': 4.229e-2/0.7, + 'X_fa': 1.546e-2/0.7, + 'X_c4': 1.561e-2/0.7, + 'X_pro': 1.511e-2/0.7, + # 'X_su': 1.1795398908390478/100, + # 'X_aa': 2.5184674870545223/100, + # 'X_fa': 0.9329873586793094/100, + # 'X_c4': 0.9313043558027272/100, + # 'X_pro': 0.3060024744735877/100, + } + +bg1 = qs.WasteStream('bg1', phase='g') +BE1 = METAB_BatchExp('BE1', outs=(bg1,), V_liq=25, V_gas=70, V_beads=10, + T=273.15+37, pH_ctrl=False, + model=adm1, fixed_headspace_P=True) +BE1.f_diff = 0.8 +BE1.set_bulk_init(**syn_ww) +BE1.set_encap_init(**h2r_seed) + +sys1 = qs.System('batch_1', path=(BE1,)) +sys1.set_dynamic_tracker(BE1, bg1) + +sys1.simulate(state_reset_hook='reset_cache', t_span=(0,8), method='BDF') +rh2 = bg1.scope.record[:, cmps.index('S_h2')] * bg1.scope.record[:, -1] # mg/L * m3/d = g(COD)/d +rh2 *= 1e3 # g(COD)/d -> mg(COD)/d +rh2 /= 1e6 # m3 reactor -> mL reactor +rh2 = rh2 * cmps.S_h2.i_mass / cmps.S_h2.chem_MW # mg/d -> mmol/d +t1 = bg1.scope.time_series +dt1 = t1[1:] - t1[:-1] +ch2 = np.cumsum(rh2[1:]*dt1) +ch2 = np.insert(ch2, 0, 0) + +#%% +tex1 = [0,1,2,3,4,5,7,9] +# yex1 = [0, 1.5541e-1, 3.2176e-2, 0, 0, 2.6512e-3, 0] # mmol/mg protein/d +yex1 = [0, 0.1554, 0.1876, 0.1226, 0.1081, 0.1108, 0.0305, 0.0738] +# yerr1 = [0, 1.2e-1, 1.33e-1, 0,0,0,0] +yerr1 = [0, 0.01231, 0.00414, 0.06947, 0.06282, 0.07872, 0.02465, 0.07875] +yex1 = [x*0.884 for x in yex1] # 88.4 mg protein/L * 10 mL = 0.884 mg +yerr1 = [x*0.884 for x in yerr1] +fig, ax = plt.subplots() +# ax.plot(t1, rh2, color='black', label='simulated') +ax.plot(t1, ch2, color='black', label='simulated') +ax.errorbar(tex1, yex1, yerr=yerr1, label='PEG', + marker='o', color='red', linestyle='dashed', + capsize=2) +ax.legend() +ax.set_xlabel('Days') +ax.set_ylabel('Cumulative H2 production [mmol]') +ax.tick_params(axis='both', which='major', direction='inout') +fig.savefig(ospath.join(figures_path, 'validation_H2E_batch.png'), + dpi=300, facecolor='white') + +#%% +# yields ~ 1.0g/L VSS +# ch4r_seed = dict( +# X_su=4.708e-1, +# X_aa=4.838e-2, +# X_fa=1.239e-1, +# X_c4=1.423e-1, +# X_pro=8.978e-2, +# X_ac=2.959e-1, +# X_h2=1.467e-1, +# # X_su=0.5, +# # X_aa=0.13, +# # X_fa=0.08, +# # X_c4=0.2, +# # X_pro=0.1, +# # X_ac=0.25, +# # X_h2=0.15 +# ) + +ch4r_seed0 = { + 'X_su': 1.956004638872764, + 'X_aa': 4.229473465645631, + 'X_fa': 1.5466419486676024, + 'X_c4': 1.5615615740437239, + 'X_pro': 0.511117652665559, + 'X_ac': 3.2327967588577007, + 'X_h2': 1.390562414197235, + # 'X_pro': 0.7, + # 'X_ac': 3.0, + # 'X_h2':1.8 + } + +T = [] +CCH4 = [] + +bg2 = qs.WasteStream('bg2', phase='g') +BE2 = METAB_BatchExp('BE2', outs=(bg2,), V_liq=25, V_gas=35, V_beads=10, + T=273.15+37, pH_ctrl=False, + # max_encapsulation_tss=6, + model=adm1, fixed_headspace_P=True) +BE2.f_diff = 0.8 +BE2.set_bulk_init(**syn_ww) + +sys2 = qs.System('batch_2', path=(BE2,)) +sys2.set_dynamic_tracker(BE2, bg2) + +# f = 0.75 +for f in [0.4, 0.75]: + ch4r_seed = {k: v*f for k,v in ch4r_seed0.items()} + BE2.set_encap_init(**ch4r_seed) + + sys2.simulate(state_reset_hook='reset_cache', t_span=(0,8), method='BDF') + # bg2.scope.plot_time_series(('S_h2', 'S_IC', 'S_ch4', 'Q')) + rch4 = bg2.scope.record[:, cmps.index('S_ch4')] * bg2.scope.record[:, -1] # mg/L * m3/d = g(COD)/d + rch4 *= 1e3 # g(COD)/d -> mg(COD)/d + rch4 /= 1e6 # m3 reactor -> mL reactor + rch4 = rch4 * cmps.S_ch4.i_mass / cmps.S_ch4.chem_MW # mg/d -> mmol/d + t2 = bg2.scope.time_series + dt2 = t2[1:] - t2[:-1] + cch4 = np.cumsum(rch4[1:]*dt2) + cch4 = np.insert(cch4, 0, 0) + T.append(t2) + CCH4.append(cch4) + + +#%% +# PEG, suspended AS +tex2 = [0,1,2,3,5,7,8] +yex2 = [0, 0.026, 0.043, 0.052, 0.067, 0.056, 0.065] +yerr2 = [0, 3.6e-3, 7.2e-3, 0.0156, 0.02, 0.013, 0.02] + +# PEG, granular AS +tex4 = [0,1,2,3,4,7,8] +yex4 = [0, 5.75e-3, 9.41e-3, 14.9e-3, 13.8e-3, 16.3e-3, 18.7e-3] +yex4 = [x*13 for x in yex4] +yerr4 = [0, 1.44e-3, 1.58e-3, 1.14e-3, 1.21e-3, 4.71e-4, 2.59e-3] +yerr4 = [x*13 for x in yerr4] + +# PEG, PAC-supported biofilm +tex5 = [0,1,2,4,6] +yex5 = [0, 1.03e-2, 2.53e-2, 8.4e-2, 9.67e-2] +yex5 = [x*17 for x in yex5] +yerr5 = [0, 3.1e-4, 3.46e-3, 6.47e-3, 7.76e-3] +yerr5 = [x*17 for x in yerr5] + +# PEG, PAC-supported biofilm, replicate +tex6 = [0,1,2,3,4,7,8] +yex6 = [0, 0.022, 0.045, 0.095, 0.121, 0.140, 0.153] +yex6 = [x*11 for x in yex6] +yerr6 = [0, 3.5e-4, 7.38e-3, 4.91e-3, 3.91e-3, 8.56e-3, 2.1e-3] +yerr6 = [x*11 for x in yerr6] + +fig, ax = plt.subplots() +# ax.plot(t2, rch4, color='black', label='simulated') +# ax.plot(t2, cch4, color='black', label='simulated') +for t, c in zip(T, CCH4): + ax.plot(t, c, color='black') +ax.errorbar(tex2, yex2, yerr=yerr2, label='PEG, suspended ADS', + marker='s', linestyle='dashed', capsize=2) +ax.errorbar(tex4, yex4, yerr=yerr4, label='PEG, granular ADS', + marker='D', linestyle='dashed', capsize=2) +ax.errorbar(tex5, yex5, yerr=yerr5, label='PEG, PAC-supported biofilm', + marker='^', linestyle='dashed', capsize=2) +ax.errorbar(tex6, yex6, yerr=yerr6, label='PEG, biofilm, replicate', + marker='o', linestyle='dashed', capsize=2) + +ax.legend(loc='upper left') +# ax.set_ylim(-0.1, 2.6) +ax.set_xlabel('Days') +ax.set_ylabel('Cumulative CH4 production [mmol]') +ax.tick_params(axis='both', which='major', direction='inout') +fig.savefig(ospath.join(figures_path, 'validation_CH4E_batch.png'), + dpi=300, facecolor='white') \ No newline at end of file diff --git a/exposan/metab_mock/__init__.py b/exposan/metab_mock/__init__.py deleted file mode 100644 index f4d991c5..00000000 --- a/exposan/metab_mock/__init__.py +++ /dev/null @@ -1,47 +0,0 @@ -# -*- coding: utf-8 -*- -''' -EXPOsan: Exposition of sanitation and resource recovery systems - -This module is developed by: - - Joy Zhang - -This module is under the University of Illinois/NCSA Open Source License. -Please refer to https://github.com/QSD-Group/EXPOsan/blob/main/LICENSE.txt -for license details. -''' - -import os - -folder = os.path.dirname(__file__) -data_path = os.path.join(folder, 'data') -results_path = os.path.join(folder, 'results') -figures_path = os.path.join(folder, 'figures') -# To save simulation results and generated figures -if not os.path.isdir(results_path): os.mkdir(results_path) -if not os.path.isdir(figures_path): os.mkdir(figures_path) - -_impact_item_loaded = False -from exposan.metab import load_lca_data - -from . import process -from .process import * - -from . import units -from .units import * - -from . import systems -from .systems import * - - -__all__ = ( - 'folder', - 'data_path', - 'results_path', - 'figures_path', - *process.__all__, - *units.__all__, - *systems.__all__, - ) - - diff --git a/exposan/metab_mock/_milestone_analysis.py b/exposan/metab_mock/_milestone_analysis.py deleted file mode 100644 index 2fea4bdb..00000000 --- a/exposan/metab_mock/_milestone_analysis.py +++ /dev/null @@ -1,126 +0,0 @@ -# -*- coding: utf-8 -*- -''' -EXPOsan: Exposition of sanitation and resource recovery systems - -This module is developed by: - - Joy Zhang - -This module is under the University of Illinois/NCSA Open Source License. -Please refer to https://github.com/QSD-Group/EXPOsan/blob/main/LICENSE.txt -for license details. -''' -from qsdsan import ImpactIndicator as IInd, ImpactItem as IItm, \ - StreamImpactItem as SIItm, \ - TEA, LCA, PowerUtility -from qsdsan.utils import ospath -from exposan.metab_mock import create_systems, data_path, results_path, DegassingMembrane -import pandas as pd - -IInd.load_from_file(ospath.join(data_path, 'TRACI_indicators.xlsx'), sheet=0) -IItm.load_from_file(ospath.join(data_path, '_impact_items.xlsx')) -IItm('Stainless_steel', source='stainless_steel') -bg_offset_CFs = IItm.get_item('biogas_offset').CFs -NaOCl_item = IItm.get_item('NaOCl') -citric_acid_item = IItm.get_item('citric_acid') - -#%% -irr = 0.1 -# p_ng = 0.65165 # USD/therm of natural gas charged to the brewery in 2016 (including tax and delivery etc.) -# p_ng = 0.85*auom('kJ').conversion_factor('therm') # [USD/kJ] 5.47 2021USD/Mcf vs. 4.19 2016USD/Mcf -op_hr = 365*24 -get = getattr -PowerUtility.price = 0.0913 # 2021$ for minnesota area - -def run_tea_lca(sys, save_report=True, info='', lt=30): - sys.simulate(state_reset_hook='reset_cache', t_span=(0, 400), method='BDF') - F = sys.flowsheet - inf = get(F.stream, f'BreweryWW_{F.ID}') - eff = get(F.stream, f'Effluent_{F.ID}') - cmps = eff.components - KJ_per_kg = cmps.i_mass/cmps.chem_MW*cmps.LHV - for ws in sys.products: - if ws.phase == 'l': - SIItm(ID=f'{ws.ID}_fugitive_ch4', linked_stream=ws, - flow_getter=lambda s: s.imass['S_ch4']*cmps.S_ch4.i_mass, - GWP100=28, MIR=0.0143794871794872) - elif ws.phase == 'g': - fget = lambda s: -sum(s.mass*KJ_per_kg)*1e-3/39 - SIItm(ID=f'{ws.ID}_NG_offset', linked_stream=ws, functional_unit='m3', - flow_getter=fget, # m3/hr natural-gas-equivalent - **bg_offset_CFs) - for u in sys.units: - if isinstance(u, DegassingMembrane): - SIItm(linked_stream=u.NaOCl, **NaOCl_item.CFs) - SIItm(linked_stream=u.citric_acid, **citric_acid_item.CFs) - - # Operation items - kWh = lambda lca: lca.system.power_utility.rate*lca.lifetime_hr - def MJ(lca): - sys = lca.system - duties = [hu.duty for hu in sys.heat_utilities] - MJ_per_hr = sum(d for d in duties if d > 0)*1e-3 - return MJ_per_hr*lca.lifetime_hr - - rCOD = inf.composite('COD', flow=True) - eff.composite('COD', flow=True) # kgCOD/hr - - tea = TEA(sys, discount_rate=irr, lifetime=lt, simulate_system=False, CEPCI=708) - levelized_cost = -tea.annualized_NPV/(rCOD * op_hr *1e-3) - print(f'{sys.ID} TEA: ${round(levelized_cost, 1)}/ton rCOD') - - lca = LCA( - sys, lifetime=lt, simulate_system=False, - electricity = kWh, - heat_onsite = MJ - ) - - gwp = lca.get_total_impacts()['GWP100'] - gwp_per_rcod = gwp/(rCOD * op_hr * lt *1e-3) - print(f'{sys.ID} LCA: {round(gwp_per_rcod, 1)} kg CO2eq/ton rCOD') - - if save_report: - tea_path = ospath.join(results_path, f'{F.ID}_{info}_tea.xlsx') - summary = { - 'lifetime':lt, - 'NPV': tea.NPV, - 'annualized_CAPEX': tea.annualized_CAPEX, - 'annualized_OPEX': tea.AOC, - 'annual utility': tea.utility_cost, - **tea.system_add_OPEX, - 'levelized_cost': levelized_cost, - 'GWP100': gwp_per_rcod - } - with pd.ExcelWriter(tea_path) as writer: - summary = pd.DataFrame.from_dict({sys.ID: summary}) - summary.to_excel(writer, sheet_name='summary') - tea.get_cashflow_table().to_excel(writer, sheet_name='cash_flow') - for u in sys.units: - rs = u.results() - rs.index.names = ['l0', 'l1'] - rs.join(pd.DataFrame.from_dict({'lifetime': u.lifetime}), - on='l1') - rs.to_excel(writer, sheet_name=u.ID) - lca.save_report(ospath.join(results_path, f'{F.ID}_{info}_lca.xlsx')) - - return sys - -#%% -if __name__ == '__main__': - # save = True - save = False - info = '35C_30yr_1yr' - lt = 30 - # sysA, sysB = create_systems(which=('A', 'B')) - # print(info) - # sysA = run_tea_lca(sysA, save, info, lt) - # run_tea_lca(sysB, save, info, lt) - systems = create_systems(which=('C',)) - for sys in systems: - u = sys.units[0] - # for tau in (1, 2): - tau = 2 - u.V_liq = tau*5 - u.V_gas = tau*5/10 - info = f'HRT{tau}' - # IItm.clear_registry() - run_tea_lca(sys, save, info, lt) diff --git a/exposan/metab_mock/data/CFs.xlsx b/exposan/metab_mock/data/CFs.xlsx deleted file mode 100644 index 5a1dcafb..00000000 Binary files a/exposan/metab_mock/data/CFs.xlsx and /dev/null differ diff --git a/exposan/metab_mock/data/HDPE_pipe_chart.xlsx b/exposan/metab_mock/data/HDPE_pipe_chart.xlsx deleted file mode 100644 index aad9e697..00000000 Binary files a/exposan/metab_mock/data/HDPE_pipe_chart.xlsx and /dev/null differ diff --git a/exposan/metab_mock/data/stainless_steel_pipe_chart.xlsx b/exposan/metab_mock/data/stainless_steel_pipe_chart.xlsx deleted file mode 100644 index 7cbc1f0c..00000000 Binary files a/exposan/metab_mock/data/stainless_steel_pipe_chart.xlsx and /dev/null differ diff --git a/exposan/metab_mock/graveyard/__init__.py b/exposan/metab_mock/graveyard/__init__.py deleted file mode 100644 index b0d9fc34..00000000 --- a/exposan/metab_mock/graveyard/__init__.py +++ /dev/null @@ -1,35 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Thu Aug 4 13:50:10 2022 - -@author: joy_c -""" - -import os -folder = os.path.dirname(__file__) -data_path = os.path.join(folder, 'data') -results_path = os.path.join(folder, 'results') -figures_path = os.path.join(folder, 'figures') -# To save simulation results and generated figures -if not os.path.isdir(results_path): os.mkdir(results_path) -if not os.path.isdir(figures_path): os.mkdir(figures_path) -del os - -from . import units -from .units import * - -from . import systems -from .systems import * - -from . import models -from .models import * - -__all__ = ( - 'folder', - 'data_path', - 'results_path', - 'figures_path', - *units.__all__, - *systems.__all__, - *models.__all__, - ) \ No newline at end of file diff --git a/exposan/metab_mock/graveyard/analyses.py b/exposan/metab_mock/graveyard/analyses.py deleted file mode 100644 index 06bdcb26..00000000 --- a/exposan/metab_mock/graveyard/analyses.py +++ /dev/null @@ -1,528 +0,0 @@ -# -*- coding: utf-8 -*- -''' -EXPOsan: Exposition of sanitation and resource recovery systems - -This module is developed by: - - Joy Zhang - -This module is under the University of Illinois/NCSA Open Source License. -Please refer to https://github.com/QSD-Group/EXPOsan/blob/main/LICENSE.txt -for license details. -''' - -from time import time -from copy import deepcopy -from qsdsan.utils import load_data, load_pickle, save_pickle, ospath -from qsdsan.stats import plot_uncertainties -from exposan.metab_mock import ( - create_systems, - create_modelA, - create_modelB, - create_modelC, - create_modelD, - create_ss_model, - run_model, - run_modelB, - run_ss_model, - results_path, - figures_path - ) -import os -import numpy as np -import pandas as pd -import matplotlib as mpl, matplotlib.pyplot as plt, matplotlib.ticker as tk - - -mpl.rcParams['font.sans-serif'] = 'arial' -mpl.rcParams["figure.autolayout"] = False -mpl.rcParams['xtick.minor.visible'] = True -mpl.rcParams['ytick.minor.visible'] = True -mpl.rcParams['savefig.facecolor'] = 'white' - -#colors -#D81B60 -#1E88E5 -#FFC107 -#004D40 - -#%% Global variables - -sysA, sysB, sysC, sysD = create_systems() -mssA = create_ss_model(sysA, R1_ID='H2E', R2_ID='CH4E') -mssC = create_ss_model(sysC, R1_ID='R1', R2_ID='R2') - -mdlA = create_modelA(sysA) -mdlB = create_modelB(sysB) -mdlC = create_modelC(sysC) -mdlD = create_modelD(sysD) - -cmps = mdlB._system.flowsheet.stream.Effluent_B.components -AnR1 = mdlB._system.flowsheet.unit.AnR1 -N = 400 -T = 200 -t_step = 5 - -keys = ['X_ch', 'S_su', 'S_va', 'S_ac', - 'X_pr', 'S_aa', 'S_bu', 'S_h2', - 'X_li', 'S_fa', 'S_pro', 'S_ch4'] -irows = [0,0,0,0,1,1,1,1,2,2,2,2] -icols = [0,1,2,3,0,1,2,3,0,1,2,3] -plots_wide = zip(keys, irows, icols) - -# indices = [(4, 5), (0, 1)] -# names = [('H2', 'CH4'), ('rCOD_1', 'rCOD_2')] -indices = [(3,4), (5,6)] -names = [('eff_COD', 'rCOD'), ('H2', 'CH4')] -pairs = zip(indices, names) - -S_keys = [ID for ID in cmps.IDs if ID.startswith('S_') and ID not in ('S_an', 'S_cat')] -X_keys = [ID for ID in cmps.IDs if ID.startswith('X_')] -gas_keys = ['S_h2_gas', 'S_ch4_gas', 'S_IC_gas'] -irows = [0,0,0,1,1,1,2,2,2,3,3,3] -icols = [0,1,2,0,1,2,0,1,2,0,1,2] -s_plots = zip(S_keys, irows, icols) -x_plots = zip(X_keys, irows, icols) -g_plots = zip(gas_keys, [0,0,0], [0,1,2]) - - -def seed_RGT(): - files = os.listdir(results_path) - seeds = [int(file_name[-3:]) for file_name in files if file_name.startswith('time_series_data')] - seed = int(str(time())[-3:]) - if len(set(seeds)) >= 1000: - raise RuntimeError('The program has run out of 3-digit seeds to use. Consider' - 'clean up the results folder.') - while seed in seeds: - seed = (seed+1) % 1000 - return seed - -#%% UA with time-series data -def single_state_var_getter(arr, unit, variable): - j = 1 - if unit in ('CH4E', 'AnR2', 'R2'): j += len(AnR1._state_keys) - if variable == 'S_IC_gas': i = 29 - elif variable == 'S_ch4_gas': i = 28 - elif variable == 'S_h2_gas': i = 27 - else: i = cmps.index(variable) - return arr[:, i+j] - -def analyze_timeseries(variable_getter, N, folder='', todf=True, **kwargs): - outputs = {} - for sample_id in range(N): - arr = np.load(ospath.join(folder, f'state_{sample_id}.npy')) - outputs[sample_id] = variable_getter(arr, **kwargs) - if todf: outputs = pd.DataFrame.from_dict(outputs) - return outputs - -def analyze_vars(seed, N, prefix='sys', sys_ID='B', pickle=True): - folder = os.path.join(results_path, f'{prefix}{sys_ID}_time_series_data_{seed}') - state_vars = S_keys+X_keys+gas_keys - dct = {} - t_arr = np.arange(0, T+t_step, t_step) - pcs = np.array([5, 25, 50, 75, 95]) - col_names = [f'{p}th percentile' for p in pcs] - if sys_ID == 'A': units = ('H2E', 'CH4E') - elif sys_ID == 'B': units = ('AnR1', 'AnR2') - else: units = ('R1', 'R2') - for u in units: - dct[u] = {} - for var in state_vars: - df = analyze_timeseries(single_state_var_getter, N=N, folder=folder, - variable=var, unit=u) - quants = np.percentile(df, q = pcs, axis=1) - df[col_names] = quants.T - df['t'] = t_arr - dct[u][var] = df - if pickle: - path = os.path.join(results_path, f'{prefix}{sys_ID}_state_vars_{seed}.pckl') - save_pickle(dct, path) - return dct - -def plot_timeseries(seed, N, unit, sys_ID='B', data=None): - if data is None: - try: - path = ospath.join(results_path, f'sys{sys_ID}_state_vars_{seed}.pckl') - data = load_pickle(path) - except: - data = analyze_vars(seed, N, sys_ID=sys_ID) - plts = deepcopy(plots_wide) - x_ticks = np.linspace(0, T, 5) - fig, axes = plt.subplots(3, 4, sharex=True, figsize=(14,8)) - for var, ir, ic in plts: - df = data[unit][var] - ax = axes[ir, ic] - # ax = axes[ic] - ys = df.iloc[:,:N].values - ax.plot(df.t, ys, color='grey', linestyle='-', alpha=0.05) - l5, l95, = ax.plot(df.t, df.loc[:,['5th percentile', '95th percentile']].values, - color='blue', linestyle='--', - label='5th, 95th') - l50, = ax.plot(df.t, df.loc[:,'50th percentile'], - color='black', linestyle='-.', - label='50th') - lct = tk.MaxNLocator(nbins=3, min_n_ticks=1) - y_ticks = lct.tick_values(np.min(ys), np.max(ys)) - ax.xaxis.set_ticks(x_ticks) - ax.yaxis.set_ticks(y_ticks) - ax.tick_params(axis='both', direction='inout', length=6, labelsize=14) - ax2x = ax.secondary_xaxis('top') - ax2x.xaxis.set_ticks(x_ticks) - ax2x.tick_params(axis='x', direction='in') - ax2x.xaxis.set_major_formatter(plt.NullFormatter()) - ax2y = ax.secondary_yaxis('right') - ax2y.yaxis.set_ticks(y_ticks) - ax2y.tick_params(axis='y', direction='in') - ax2y.yaxis.set_major_formatter(plt.NullFormatter()) - # if ir == 0 and ic == 0: ax.legend(handles=[l5, l50, lbl]) - plt.ylim(y_ticks[0], y_ticks[-1]) - plt.subplots_adjust(wspace=0.2, hspace=0.1) - fig.savefig(ospath.join(figures_path, f'{unit}_states.png'), dpi=300) - return fig, axes - -def plot_kde2d_metrics(seed, model, sys_ID): - metrics = model.metrics - if model.table is None: - path = ospath.join(results_path, f'sys{sys_ID}_table_{seed}.xlsx') - model.table = load_data(path=path, header=[0,1]) - for ids, names in pairs: - i, j = ids - x, y = names - fig, ax = plot_uncertainties(model, x_axis=metrics[i], y_axis=metrics[j], - kind='kde-box', center_kws={'fill':True}, - margin_kws={'width':0.5}) - ax0, ax1, ax2 = fig.axes # KDE, top box, right box - ax0x = ax0.secondary_xaxis('top') - ax0y = ax0.secondary_yaxis('right') - for ax in (ax0, ax0x, ax0y): - ax.tick_params(axis='both', which='both', direction='inout', width=1) - - for txt in ax0.get_xticklabels()+ax0.get_yticklabels(): - txt.set_fontsize(16) - ax0x.set_xticklabels('') - ax0y.set_yticklabels('') - - ax1.xaxis.set_visible(False) - ax1.spines.clear() - ax2.spines.clear() - ax2.yaxis.set_visible(False) - xl, xu = ax0.get_xlim() - yl, yu = ax0.get_ylim() - if yl < 0: - ax0.set_ylim((0, yu)) - yl, yu = ax0.get_ylim() - ax0.set_xlim((xl, xu)) - ax0.set_ylim((yl, yu)) - fig.savefig(ospath.join(figures_path, f'sys{sys_ID}_{x}_vs_{y}.png'), dpi=300) - del fig, ax - -def plot_scatter(seed, model1, model2=None, bl_metrics=None): - if model2 is None: - mdls = (model1,) - mks = ('o',) - clrs = ('#E66100',) - prefix = model1.system.flowsheet.ID[-1] - else: - mdls = (model1, model2) - mks = ('o', '^') - clrs = ('#D81B60', '#1E88E5') - prefix = f'{model1.system.flowsheet.ID[-1]}v{model2.system.flowsheet.ID[-1]}' - nx = len(model1.parameters) - ny = len(model1.metrics) - for mdl in mdls: - if mdl.table is None: - path = ospath.join(results_path, f'sys{mdl.system.flowsheet.ID[-1]}_table_{seed}.xlsx') - mdl.table = load_data(path=path, header=[0,1]) - dfs_x = [mdl.table.iloc[:, :nx] for mdl in mdls] - dfs_y = [mdl.table.iloc[:, nx:] for mdl in mdls] - - fig, axes = plt.subplots(ny, nx, sharex=False, sharey=False, - figsize=(nx*2.5, ny*2.5), - layout='constrained') - xlct = tk.MaxNLocator(nbins=3, min_n_ticks=1) - ylct = tk.MaxNLocator(nbins=3, min_n_ticks=1) - for j in range(ny): - ys = [df_y.iloc[:,j].values for df_y in dfs_y] - try: yb = bl_metrics[j] - except: yb = -1 - ymin = max(0, min(np.min(ys), yb)) - ymax = max(np.max(ys), yb) - y_ticks = ylct.tick_values(ymin, ymax) - for i in range(nx): - xs = [df_x.iloc[:,i].values for df_x in dfs_x] - x_ticks = xlct.tick_values(np.min(xs), np.max(xs)) - if nx > 1: ax = axes[j,i] - else: ax = axes[j] - ax.axhline(y=yb, ls='-', lw=0.5, c='black') - for x, y, m, c in zip(xs, ys, mks, clrs): - ax.scatter(x, y, marker=m, s=1, c=c) - ax.tick_params(axis='both', which='both', - direction='inout', labelsize=10.5) - ax.xaxis.set_ticks(x_ticks) - ax.yaxis.set_ticks(y_ticks) - if i > 0: ax.yaxis.set_ticklabels([]) - if j < ny-1: ax.xaxis.set_ticklabels([]) - ax2x = ax.secondary_xaxis('top') - ax2x.xaxis.set_ticks(x_ticks) - ax2x.tick_params(axis='x', which='both', direction='in') - ax2x.xaxis.set_ticklabels([]) - ax2y = ax.secondary_yaxis('right') - ax2y.yaxis.set_ticks(y_ticks) - ax2y.tick_params(axis='y', which='both', direction='in') - ax2y.yaxis.set_ticklabels([]) - fig.savefig(ospath.join(figures_path, f'{prefix}_table_{seed}.png'), dpi=300) - return fig, axes - -def plot_2d_scatter(seed, model, bl_metrics=None): - npr = len(model.parameters) - nmt = len(model.metrics) - if model.table is None: - path = ospath.join(results_path, f'sysC_table_{seed}.xlsx') - model.table = load_data(path=path, header=[0,1]) - df_p = model.table.iloc[:, :npr] - df_m = model.table.iloc[:, npr:] - nplot = ceil(nmt/2) - nrow = ceil(nplot/2) - xlct = tk.MaxNLocator(nbins=3, min_n_ticks=1) - ylct = tk.MaxNLocator(nbins=3, min_n_ticks=1) - for i in range(npr): - fig, axes = plt.subplots(nrow, 2, sharex=False, sharey=False, - figsize=(6, nrow*2.5)) - c = df_p.iloc[:,i].values - for j, ax in zip(range(nplot), axes.ravel()): - x = df_m.iloc[:,j*2].values - y = df_m.iloc[:,j*2+1].values - try: xmin = np.min(x) - except: breakpoint() - xmax = np.max(x) - ymin = np.min(y) - ymax = np.max(y) - im = ax.scatter(x, y, marker='o', s=1, c=c, cmap='viridis') - if bl_metrics is not None: - x_bl, y_bl = bl_metrics[j*2: (j*2+2)] - ax.scatter(x_bl, y_bl, s=1.5, marker='^', c='black') - if x_bl < xmin: xmin = x_bl - elif x_bl > xmax: xmax = x_bl - if y_bl < ymin: ymin = y_bl - elif y_bl > ymax: ymax = y_bl - ax.tick_params(axis='both', which='both', - direction='inout', labelsize=10.5) - x_ticks = xlct.tick_values(xmin, xmax) - y_ticks = ylct.tick_values(ymin, ymax) - ax.xaxis.set_ticks(x_ticks) - ax.yaxis.set_ticks(y_ticks) - ax2x = ax.secondary_xaxis('top') - ax2x.xaxis.set_ticks(x_ticks) - ax2x.tick_params(axis='x', which='both', direction='in') - ax2x.xaxis.set_ticklabels([]) - ax2y = ax.secondary_yaxis('right') - ax2y.yaxis.set_ticks(y_ticks) - ax2y.tick_params(axis='y', which='both', direction='in') - ax2y.yaxis.set_ticklabels([]) - cbar = fig.colorbar(im, ax=ax) - cbar.ax.tick_params(labelsize=10.5) - fig.subplots_adjust(wspace=0., hspace=0.) - fig.savefig(ospath.join(figures_path, f'2dscatter_param{i}_{seed}.png'), dpi=300) - del fig, axes - -def plot_ss_convergence(seed, N, unit='CH4E', prefix='ss', sys_ID='A', data=None, - baseline_sys=None, baseline_unit_ID='AnR2'): - if data is None: - try: - path = ospath.join(results_path, f'{prefix}{sys_ID}_state_vars_{seed}.pckl') - data = load_pickle(path) - except: - data = analyze_vars(seed, N, sys_ID=sys_ID) - if baseline_sys is not None: - blu = getattr(baseline_sys.flowsheet.unit, baseline_unit_ID) - bl_df = pd.DataFrame(blu.scope.record, columns=blu._state_keys) - bl_ts = blu.scope.time_series - fig, axes = plt.subplots(4, 3, sharex=True, figsize=(12,12)) - x_ticks = np.linspace(0, 200, 6) - for group, par_size in [(s_plots, 's'), (x_plots, 'x'), (g_plots, 'g')]: - if par_size == 'g': - fig, axes = plt.subplots(1, 3, sharex=True, figsize=(12,3.2), layout='constrained') - else: - fig, axes = plt.subplots(4, 3, sharex=True, figsize=(12,12), layout='constrained') - plts = deepcopy(group) - for var, ir, ic in plts: - df = data[unit][var] - if par_size == 'g': ax = axes[ic] - else: ax = axes[ir, ic] - ys = df.iloc[:,:N].values - ax.plot(df.t, ys, color='grey', linestyle='-', alpha=0.25) - if baseline_sys is not None: - bl_y = bl_df.loc[:, var].values - ax.plot(bl_ts, bl_y, color='#FFC107') - ymin = np.min(bl_y) - ymax = np.max(bl_y) - else: - ymin = np.inf - ymax = -1 - lct = tk.MaxNLocator(nbins=3, min_n_ticks=1) - y_ticks = lct.tick_values(min(np.min(ys), ymin), max(np.max(ys), ymax)) - ax.xaxis.set_ticks(x_ticks) - ax.yaxis.set_ticks(y_ticks) - ax.tick_params(axis='both', which='both', - direction='inout', labelsize=14) - ax2x = ax.secondary_xaxis('top') - ax2x.xaxis.set_ticks(x_ticks) - ax2x.tick_params(axis='x', which='both', direction='in') - ax2x.xaxis.set_ticklabels([]) - ax2y = ax.secondary_yaxis('right') - ax2y.yaxis.set_ticks(y_ticks) - ax2y.tick_params(axis='y', which='both', direction='in') - ax2y.yaxis.set_ticklabels([]) - del plts - # fig.subplots_adjust(hspace=0., wspace=0.) - fig.savefig(ospath.join(figures_path, f'{prefix}{sys_ID}_{par_size}t_{unit}.png'), dpi=300) - # fig.savefig(ospath.join(figures_path, f'{prefix}{sys_ID}_{par_size}t_tau01.png'), dpi=300) - del fig, axes - -def plot_area(seed, model): - if model.table is None: - path = ospath.join(results_path, f'sys{model.system.flowsheet.ID[-1]}_table_{seed}.xlsx') - model.table = load_data(path=path, header=[0,1]) - data = model.table - data.columns = ['x', 'r1', 'rt', 'h2cod', 'ch4cod'] - data.sort_values(by='x', axis=0, inplace=True) - x = data.iloc[:,0].to_numpy(copy=True) - y_r1 = data.iloc[:,1].to_numpy(copy=True) - y_tot = data.iloc[:,2].to_numpy(copy=True) - fig, ax = plt.subplots(figsize=(4,3), layout='constrained') - ax.plot(x, y_r1, '-', c='#D81B60', linewidth=0.4) - area1 = ax.fill_between(x, 0., y_r1, color='#D81B60', alpha=0.7, label='Stage 1') - ax.plot(x, y_tot, '-', c='#1E88E5', linewidth=0.4) - area2 = ax.fill_between(x, y_r1, y_tot, color='#1E88E5', alpha=0.7, label='Stage 2') - ax.set_xlim(np.min(x), np.max(x)) - ax.set_ylim(bottom=0., top=100) - # ax.set_yscale('log') - ax.set_xlabel('Recirculation rate [-]', fontsize=11) - ax.set_ylabel('COD removal [%]', fontsize=11) - ax.legend(handles=[area1, area2], fontsize=10) - ax.tick_params(axis='both', which='both', direction='inout', labelsize=10) - ax2x = ax.secondary_xaxis('top') - ax2x.set_xlim(np.min(x), np.max(x)) - ax2x.tick_params(direction='in', which='both') - ax2x.xaxis.set_ticklabels([]) - ax2y = ax.secondary_yaxis('right') - ax2y.set_ylim(bottom=0., top=100) - # ax2y.set_yscale('log') - ax2y.tick_params(direction='in', which='both') - ax2y.yaxis.set_ticklabels([]) - fig.savefig(ospath.join(figures_path, f'rCOD_vs_recir_{seed}.png'), dpi=300) - return fig, ax - -#%% -import seaborn as sns -from math import ceil - -def f_grouping(cod_removal): - if cod_removal > 50: return 'high' - else: return 'low' - -def plot_dist_bygroup(seed, model, sys_ID='A'): - group = model.table.loc[:,('System','Total COD removal [%]')].apply(f_grouping) - n_param = len(model.parameters) - ncol = 6 - nrow = ceil(n_param/ncol) - x_df = model.table.iloc[:, :n_param] - fig, axes = plt.subplots(nrow, ncol, sharey=True, figsize=(ncol*3, nrow*2.5)) - for col, ax in zip(x_df, axes.ravel()): - unit = col[0].split('-')[-1] - var = col[1].split(' 0')[0] - sns.histplot(data=x_df, x=col, hue=group, bins=30, stat='probability', - cumulative=True, - common_norm=False, kde=True, legend=False, ax=ax) - ax.set_xlabel(f'{unit} {var}_0') - fig.subplots_adjust(hspace=0., wspace=0.) - fig.savefig(ospath.join(figures_path, f'hist_ss{sys_ID}_{seed}'), dpi=300) - - -#%% -def run_UA_AvC(seed=None, N=N, T=T, t_step=t_step, plot=True): - seed = seed or seed_RGT() - run_model(mdlA, N, T, t_step, method='BDF', seed=seed) - run_model(mdlC, N, T, t_step, method='BDF', seed=seed) - # outC = analyze_vars(seed, N, prefix='sys', sys_ID='C') - for p in mdlB.parameters: - p.setter(p.baseline) - temp_sysB = mdlB._system - temp_sysB.simulate(t_span=(0, T), method='BDF', - state_reset_hook='reset_cache') - tmB = [m.getter() for m in mdlB.metrics] - if plot: - plot_scatter(seed, mdlA, mdlC, tmB) - # plot_ss_convergence(seed, N, unit='R1', prefix='sys', sys_ID='C', data=outC, - # baseline_sys=temp_sysB, baseline_unit_ID='AnR1') - # plot_ss_convergence(seed, N, unit='R2', prefix='sys', sys_ID='C', data=outC, - # baseline_sys=temp_sysB, baseline_unit_ID='AnR2') - print(f'Seed used for uncertainty analysis of systems A & C is {seed}.') - for mdl in (mdlA, mdlC): - for p in mdl.parameters: - p.setter(p.baseline) - return seed - -def run_UA_DvB(seed=None, N=N, T=T, t_step=t_step, run=True, plot=True): - seed = seed or seed_RGT() - if run: run_model(mdlD, N, T, t_step, method='BDF', seed=seed) - for p in mdlB.parameters: - p.setter(p.baseline) - temp_sysB = mdlB._system - temp_sysB.simulate(t_span=(0, T), method='BDF', - state_reset_hook='reset_cache') - tmB = [m.getter() for m in mdlB.metrics] - if plot: - plot_scatter(seed, mdlD, None, tmB) - # plot_ss_convergence(seed, N, unit='R1', prefix='sys', sys_ID='D', - # baseline_sys=temp_sysB, baseline_unit_ID='AnR1') - print(f'Seed used for uncertainty analysis of system D is {seed}.') - for p in mdlD.parameters: - p.setter(p.baseline) - return seed - -def run_UA_sysB(seed=None, N=N, T=T, t_step=t_step, plot=True): - seed = seed or seed_RGT() - run_modelB(mdlB, N, T, t_step, method='BDF', seed=seed) - out = analyze_vars(seed, N, sys='B') - if plot: - for u in ('AnR1', 'AnR2'): - plot_timeseries(seed, N, u, data=out, sys='B') - plot_kde2d_metrics(seed, mdlB, sys_ID='B') - print(f'Seed used for uncertainty analysis of system B is {seed}.') - for p in mdlB.parameters: - p.setter(p.baseline) - return seed - -def run_ss_AvC(seed=None, N=N, T=T, t_step=t_step, plot=True): - seed = seed or seed_RGT() - run_ss_model(mssA, N, T, t_step, sys_ID='A', R1_ID='H2E', R2_ID='CH4E', seed=seed) - run_ss_model(mssC, N, T, t_step, sys_ID='C', R1_ID='R1', R2_ID='R2', seed=seed) - outA = analyze_vars(seed, N, prefix='ss', sys_ID='A') - outC = analyze_vars(seed, N, prefix='ss', sys_ID='C') - if plot: - plot_ss_convergence(seed, N, unit='CH4E', prefix='ss', sys_ID='A', data=outA) - plot_ss_convergence(seed, N, unit='R2', prefix='ss', sys_ID='C', data=outC) - plot_dist_bygroup(seed, mssA) - plot_dist_bygroup(seed, mssC, sys_ID='C') - print(f'Seed used for convergence test is {seed}.') - return seed - -def dv1_analysis(seed=None, N=100, T=T, t_step=t_step, run=True, plot=True): - seed = seed or seed_RGT() - if run: run_model(mdlD, N, T, t_step, method='BDF', seed=seed) - if plot: fig, ax = plot_area(seed, mdlD) - return seed - -#%% -if __name__ == '__main__': - # run_UA_AvC(seed=628) - # run_UA_DvB(seed=96) - # run_ss_AvC(seed=223, N=100) - dv1_analysis(seed=796) - # seed = 952 - # run_modelB(mdl, N, T, t_step, method='BDF', seed=seed) - # out = analyze_vars(seed, N, 'B') - # plot_timeseries(seed, N, 'AnR1', 'B') - # plot_timeseries(seed, N, 'AnR2', 'B') - # plot_kde2d_metrics(seed, 'B') diff --git a/exposan/metab_mock/graveyard/degas_test.py b/exposan/metab_mock/graveyard/degas_test.py deleted file mode 100644 index 0340b4d1..00000000 --- a/exposan/metab_mock/graveyard/degas_test.py +++ /dev/null @@ -1,344 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Tue Oct 25 10:39:19 2022 - -@author: joy_c -""" -# from qsdsan import Components, set_thermo, Process, Processes, WasteStream, System -import qsdsan as qs, qsdsan.processes as pc, qsdsan.sanunits as su, numpy as np, \ - matplotlib.pyplot as plt, matplotlib as mpl#, matplotlib.ticker as tk -from qsdsan.utils import ospath, load_data -from exposan.metab_mock import DegassingMembrane as DM, results_path, figures_path -from chaospy import distributions as shape - -def create_cmps(): - adm_cmps = pc.create_adm1_cmps(False) - cmps = qs.Components([adm_cmps.S_h2, adm_cmps.S_ch4, adm_cmps.S_IC, adm_cmps.H2O]) - cmps.default_compile() - qs.set_thermo(cmps) - return cmps - -R = 8.3145e-2 # Universal gas constant, [bar/M/K] -T_op = 298.15 -rhos = np.zeros(6) -rhos[:3] = [1e-3,1e-2,2e-2] - -def f_rhos(state_arr, params): - # S_h2, S_ch4, S_IC, H2O, g_h2, g_ch4, g_co2, Q = state_arr - KH = params['K_H'] - kLa = params['kLa'] - biogas_S = state_arr[:3].copy() - biogas_p = R * T_op * state_arr[4:7] - rhos[-3:] = kLa * (biogas_S - KH * biogas_p) - rhos[0] = params['rh2'] - return rhos - -def create_processes(): - p1 = qs.Process('H2_production', - reaction={'S_h2':1}, - ref_component='S_h2', - conserved_for=()) - p2 = qs.Process('CH4_production', - reaction={'S_ch4':1}, - ref_component='S_ch4', - conserved_for=()) - p3 = qs.Process('CO2_production', - reaction={'S_IC':1}, - ref_component='S_IC', - conserved_for=()) - p4 = qs.Process('H2_transfer', - reaction={'S_h2':-1}, - ref_component='S_h2', - conserved_for=()) - p5 = qs.Process('CH4_transfer', - reaction={'S_ch4':-1}, - ref_component='S_ch4', - conserved_for=()) - p6 = qs.Process('CO2_transfer', - reaction={'S_IC':-1}, - ref_component='S_IC', - conserved_for=()) - - ps = qs.Processes([p1,p2,p3,p4,p5,p6]) - ps.compile() - ps.set_rate_function(f_rhos) - ps.rate_function._params = { - 'K_H': np.array([7.8e-4, 1.4e-3, 3.5e-2]), - 'kLa':200, - 'rh2':1e-3} - - ps.__dict__['_biogas_IDs'] = ('S_h2', 'S_ch4', 'S_IC') - return ps - -Q = 5 # influent flowrate [m3/d] -T1 = 273.15+35 # temperature [K] -Vl1 = 5 # liquid volume [m^3] -Vg1 = 0.556 # headspace volume [m^3] -tau_1 = 0.021 # degassing membrane retention time [d] -f_Q = 1 # recirculation rate - -def create_system(fixed_hsp_P=False): - cmps = create_cmps() - ps = create_processes() - inf = qs.WasteStream('inf') - eff = qs.WasteStream('eff') - bgh = qs.WasteStream('bgh', phase='g') - bgm = qs.WasteStream('bgm', phase='g') - inf.set_flow_by_concentration(Q, {}, units=('m3/d', 'kg/m3')) - R1 = su.AnaerobicCSTR('R1', ins=[inf, 'return_1'], - outs=(bgh, 'sidestream_1', eff), - T=298.15, V_liq=Vl1, V_gas=Vg1, model=ps, - split=(f_Q, 1), headspace_P=1.013, - fixed_headspace_P=fixed_hsp_P) - R1.set_init_conc(S_h2=8.503, S_ch4=0.0422, S_IC=1141.2) - DM1 = DM('DM1', ins=R1-1, outs=(bgm, 1-R1), tau=tau_1) - sys = qs.System('degas', path=(R1, DM1), recycle=(DM1-1,)) - return sys - -#%% -def add_params(model, fix_H2_rprod=False): - param = model.parameter - R1 = model._system.flowsheet.unit.R1 - ps = R1.model - - b = 1 - D = shape.Uniform(1, 9) - @param(name='Recirculation_rate', element=R1, kind='coupled', - units='', baseline=b, distribution=D) - def set_fQ(fq): - R1.split = [fq, 1] - - if not fix_H2_rprod: - b = 1e-3 - D = shape.LogUniform(-6*np.log(10), -0.1*np.log(10)) - @param(name='H2_production_rate', element=ps, kind='coupled', - units='kgCOD/m3/d', baseline=b, distribution=D) - def set_H2_rprod(r): - ps.rate_function._params['rh2'] = r - - return model - -KI_h2_fa=5e-6 -KI_h2_c4=1e-5 -KI_h2_pro=3.5e-6 -def get_I(R1, Ki): - Si = R1._state[0] - return pc.non_compet_inhibit(Si, Ki) - -def add_metrics(model): - metric = model.metric - sreg = model._system.flowsheet.stream - bgm, bgh, eff = sreg.bgm, sreg.bgh, sreg.eff - R1 = model._system.flowsheet.unit.R1 - h2_i_mass = eff.components.S_h2.i_mass - - @metric(name='Sidestream_H2_extraction', units='kg/d', element='Biogas') - def get_ss_H2(): - return bgm.imass['S_h2'] * h2_i_mass * 24 - - @metric(name='Headspace_H2_extraction', units='kg/d', element='Biogas') - def get_hs_H2(): - return bgh.imass['S_h2'] * h2_i_mass * 24 - - @metric(name='Dissolved_H2', units='kgCOD/m3', element='Reactor') - def get_S_h2(): - return eff.iconc['S_h2']/1e3 - - @metric(name='H2_inhibition_fa', units='', element='Reactor') - def get_Ifa(): - return get_I(R1, KI_h2_fa) - - @metric(name='H2_inhibition_c4', units='', element='Reactor') - def get_Ic4(): - return get_I(R1, KI_h2_c4) - - @metric(name='H2_inhibition_pro', units='', element='Reactor') - def get_Ipro(): - return get_I(R1, KI_h2_pro) - - return model - - - -#%% -mpl.rcParams['font.sans-serif'] = 'arial' -mpl.rcParams["figure.autolayout"] = False -mpl.rcParams['xtick.minor.visible'] = True -mpl.rcParams['ytick.minor.visible'] = True -mpl.rcParams['savefig.facecolor'] = 'white' - -def meshgrid_sample(p1, p2, n): - xl, xu = p1.distribution.lower[0], p1.distribution.upper[0] - yl, yu = np.log10(p2.distribution.lower[0]), np.log10(p2.distribution.upper[0]) - x = np.linspace(xl, xu, n) - y = 10**np.linspace(yl, yu, n) - xx, yy = np.meshgrid(x, y) - samples = np.array([xx.reshape(n**2), yy.reshape(n**2)]) - return samples.T, xx, yy - -def run_mapping(model, n, T, t_step, run=True, method='BDF', mpath=''): - x = model.parameters[0] - y = model.parameters[1] - samples, xx, yy = meshgrid_sample(x, y, n) - if run: - model.load_samples(samples) - t_span = (0, T) - t_eval = np.arange(0, T+t_step, t_step) - model.evaluate( - state_reset_hook='reset_cache', - t_span=t_span, - t_eval=t_eval, - method=method - ) - if mpath: model.table.to_excel(mpath) - return xx, yy - -def fmt(x): - s = f"{x:.1e}" - return rf"{s}" if plt.rcParams["text.usetex"] else f"{s}" - -irs = [0, 0, 1, 1, 2, 2] -ics = [0, 1, 0, 1, 0, 1] - -def plot_heatmaps(xx, yy, model=None, path='', wide=False): - if model: data = model.table - else: - path = path or ospath.join(results_path, 'table_2dv.xlsx') - data = load_data(path, sheet=0, header=[0,1]) - zs = data.iloc[:,2:].to_numpy(copy=True) - n = int(zs.shape[0] ** 0.5) - zs = zs.T - zz = zs.reshape((zs.shape[0], n, n)) - if wide: - plts = zip(zz, ics, irs) - fig, axes = plt.subplots(2, 3, sharex=True, sharey=True, figsize=(14,9)) - else: - plts = zip(zz, irs, ics) - fig, axes = plt.subplots(3, 2, sharex=True, sharey=True, figsize=(9,12)) - - fig.set_tight_layout({'h_pad':4, 'w_pad': 0}) - for z, ir, ic in plts: - ax = axes[ir, ic] - ax.set_yscale('log') - if ic > 0: - norm = 'linear' - cfmt = lambda x: f"{x:.2g}" - else: - norm = 'log' - cfmt = fmt - pos = ax.pcolormesh(xx, yy, z, norm=norm, shading='gouraud') - cbar = fig.colorbar(pos, ax=ax) - cbar.ax.tick_params(labelsize=14) - ax.tick_params(axis='both', which='both', direction='inout', labelsize=14) - ax2x = ax.secondary_xaxis('top') - ax2x.tick_params(direction='in', which='both') - ax2x.xaxis.set_ticklabels([]) - ax2y = ax.secondary_yaxis('right') - ax2y.tick_params(direction='in', which='both') - ax2y.yaxis.set_ticklabels([]) - cs = ax.contour(xx, yy, z, - colors='white', norm=norm, origin='lower', - linestyles='dashed', linewidths=1, - extent=(xx[0,0], xx[0,-1], yy[0,0], yy[-1,0])) - ax.clabel(cs, cs.levels, inline=True, fmt=cfmt, fontsize=11) - fig.savefig(ospath.join(figures_path, 'heatmaps.png'), dpi=300) - - -def dv_analysis(n=20, T=100, t_step=10, run=True, save_to='table_2dv.xlsx', - plot=True, wide=True): - path = ospath.join(results_path, save_to) - sys = create_system() - mdl = qs.Model(sys) - mdl = add_params(mdl) - mdl = add_metrics(mdl) - xx, yy = run_mapping(mdl, n, T, t_step, run, mpath=path) - if plot: - if run: - plot_heatmaps(xx, yy, mdl, wide=wide) - for p in mdl.parameters: - p.setter(p.baseline) - else: - plot_heatmaps(xx, yy, path=path, wide=wide) - return xx, yy - -#%% -import pandas as pd - -def run_scenarios(model, N=50, - T=100, t_step=10, method='BDF', mpath=''): - X, = model.parameters - xl, = X.distribution.lower - xu, = X.distribution.upper - samples = np.array([np.linspace(xl, xu, N),]).T - model.load_samples(samples) - t_span = (0, T) - t_eval = np.arange(0, T+t_step, t_step) - mpath = mpath or ospath.join(results_path, 'table_1dv.xlsx') - ps = model._system.flowsheet.unit.R1.model - dct = {} - for rh2 in 10**np.linspace(-6,0,7): - ps.rate_function._params['rh2'] = rh2 - model.evaluate( - state_reset_hook='reset_cache', - t_span=t_span, - t_eval=t_eval, - method=method - ) - plot_area(model.table, rh2) - dct[f'{rh2:.0e}'] = model.table.copy() - writer=pd.ExcelWriter(mpath) - for k, df in dct.items(): - df.to_excel(writer, sheet_name=k) - writer.save() - -def plot_area(data, rh2): - name = f'{rh2:.0e}' - x = data.iloc[:,0].to_numpy(copy=True) - y_mem = data.iloc[:,1].to_numpy(copy=True) - y_tot = y_mem + data.iloc[:,2].to_numpy() - fig, ax = plt.subplots(figsize=(4,3), layout='constrained') - ax.plot(x, y_mem, '-', c='#D81B60', linewidth=0.4) - area1 = ax.fill_between(x, 0., y_mem, color='#D81B60', alpha=0.7, label='Sidestream') - ax.plot(x, y_tot, '-', c='#1E88E5', linewidth=0.4) - area2 = ax.fill_between(x, y_mem, y_tot, color='#1E88E5', alpha=0.7, label='Headspace') - ax.set_ylim(bottom=0., top=rh2*0.7) - ax.set_xlim(np.min(x), np.max(x)) - ax.set_xlabel('Recirculation rate [-]', fontsize=11) - ax.set_ylabel('H2 extraction rate [kg/d]', fontsize=11) - ax.legend(handles=[area1, area2], fontsize=10) - ax.set_title(f'H2 production rate = {name} kgCOD/m3/d', fontsize=11) - ax.tick_params(axis='both', which='both', direction='inout', labelsize=10) - ax2x = ax.secondary_xaxis('top') - ax2x.set_xlim(np.min(x), np.max(x)) - ax2x.tick_params(direction='in', which='both') - ax2x.xaxis.set_ticklabels([]) - ax2y = ax.secondary_yaxis('right') - ax2y.set_ylim(bottom=0., top=rh2*0.7) - ax2y.tick_params(direction='in', which='both') - ax2y.yaxis.set_ticklabels([]) - fig.savefig(ospath.join(figures_path, f'area_rH2_{name}.png'), dpi=300) - del fig, ax - -def dv1_analysis(N=50, T=100, t_step=10, run=True, fixed_hsp_P=False, - save_to='table_1dv_fixedP.xlsx', plot=True): - path = ospath.join(results_path, save_to) - if run: - sys = create_system(fixed_hsp_P) - mdl2 = qs.Model(sys) - mdl2 = add_params(mdl2, True) - mdl2 = add_metrics(mdl2) - run_scenarios(mdl2, N=N, T=T, t_step=t_step, mpath=path) - elif plot: - dct = load_data(path, sheet=None, header=[0,1], - skiprows=[2,], index_col=0) - for k, v in dct.items(): - rh2 = float(k) - plot_area(v, rh2) - - -#%% -if __name__ == '__main__': - # xx, yy = dv_analysis(n=4, run=False) - # xx, yy = dv_analysis(run=False) - # xx, yy = dv_analysis() - dv1_analysis(save_to='table_1dv_Pvar.xlsx') \ No newline at end of file diff --git a/exposan/metab_mock/graveyard/individual_test.py b/exposan/metab_mock/graveyard/individual_test.py deleted file mode 100644 index f26ce779..00000000 --- a/exposan/metab_mock/graveyard/individual_test.py +++ /dev/null @@ -1,333 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Oct 21 21:19:00 2022 - -@author: yalinli_cabbi -""" -import numpy as np, pandas as pd, qsdsan as qs -from exposan.metab_mock import ( - create_modelB, - create_modelC, - create_systems, - R1_ss_conds, - R2_ss_conds, - ) - -mdlB = create_modelB() -sysB = mdlB.system -sysB.simulate(state_reset_hook='reset_cache', t_span=(0,200), method='BDF') - - -mdlC = create_modelC() -sysC = mdlC.system -fC = sysC.flowsheet -inf = fC.stream.BreweryWW_C -eff = fC.stream.Effluent_C -bgs = [s for s in sysC.products if 'biogas' in s.ID] - -R1 = fC.unit.R1 -R2 = fC.unit.R2 -R1.set_init_conc(**R1_ss_conds) -R2.set_init_conc(**R2_ss_conds) - -DM1_c = fC.unit.DM1_c -DM2_c = fC.unit.DM2_c - -for u in (DM1_c, DM2_c): - for gas in ('H2', 'CH4', 'CO2'): - setattr(u, f'{gas}_degas_efficiency', 1) - -cmps = qs.get_components() -def get_COD(streams): - try: iter(streams) - except: streams = [streams] - return sum((s.mass*cmps.i_COD).sum() for s in streams) - -# # Parameter names -# H2E sidestream split -# CH4E sidestream split -# H2 removal efficiency -# CH4 removal efficiency -# CO2 removal efficiency - -def get_p(mdl, name): - for p in mdlC.parameters: - if p.name == name: return p -p_R1_split = get_p(mdlC, 'H2E sidestream split') -p_R2_split = get_p(mdlC, 'CH4E sidestream split') - -# # Metric names -# SRT -# R1 VFAs -# R1 H2 -# R2 CH4 -# Effluent COD -# Total COD removal -# H2 production -# CH4 production - -def get_m(mdl, name): - for m in mdl.metrics: - if m.name == name: return m -mB = get_m(mdlB, 'Total COD removal') -mC = get_m(mdlC, 'Total COD removal') - - -dcts = cod_rm_dct, cod_eff_dct, cod_bgs_dct = [{}, {}, {}] -splits = [round(i, 1) for i in np.arange(0.1, 1, 0.1)] -for R1_split in splits: - p_R1_split.setter(R1_split) - cod_rm_dct[R1_split] = [] - cod_eff_dct[R1_split] = [] - cod_bgs_dct[R1_split] = [] - for R2_split in splits: - p_R2_split.setter(R2_split) - R1.set_init_conc(**R1_ss_conds) - R2.set_init_conc(**R2_ss_conds) - try: - sysC.simulate(state_reset_hook='reset_cache', t_span=(0,200), method='BDF') - cod_rm_dct[R1_split].append(mC()) - cod_eff_dct[R1_split].append(get_COD(eff)) - cod_bgs_dct[R1_split].append(get_COD(bgs)) - except: - for dct in dcts: dct[R1_split].append(np.nan) - -dfs = [pd.DataFrame.from_dict(dct) for dct in dcts] -for df in dfs: - df.columns.name = 'R1 split' - df.index = splits - df.index.name = 'R2 split' - - -# R1_split = 0.1 -# p_R1_split.setter(R1_split) -# lst = [] -# splits = [round(i, 1) for i in np.arange(0.1, 1, 0.1)] -# for R2_split in splits: -# p_R2_split.setter(R2_split) -# R1.set_init_conc(**R1_ss_conds) -# R2.set_init_conc(**R2_ss_conds) -# try: -# sysC.simulate(state_reset_hook='reset_cache', t_span=(0,200), method='BDF') -# removal = m() -# except: removal = np.nan -# print(removal) -# lst.append(removal) - - -# %% - -import numpy as np, pandas as pd, qsdsan as qs -from exposan.metab_mock import create_systems, create_modelC - -sysA, sysB, sysC = create_systems() - -fB = sysB.flowsheet -fC = sysC.flowsheet -inf = fC.stream.BreweryWW_C -cmps = inf.components -eff = fC.stream.Effluent_C -bgs = [s for s in sysC.products if 'biogas' in s.ID] - -def get_COD(streams): - try: iter(streams) - except: streams = [streams] - return sum((s.mass*cmps.i_COD).sum() for s in streams) - -R1 = fC.unit.R1 -R2 = fC.unit.R2 - -R2.split = (0.1, 0.9) - -splits = np.linspace(0.1, 0.9, 9) - -inf_CODs = [] -eff_CODs = [] -bgs_CODs = [] -for split in splits: - R1.split = (split, 1-split) - # sysC.set_tolerance(rmol=1e-6) - try: - sysC.simulate(state_reset_hook='reset_cache', t_span=(0,200), method='BDF') - inf_CODs.append(get_COD(inf)) - eff_CODs.append(get_COD(eff)) - bgs_CODs.append(get_COD(bgs)) - except: - inf_CODs.append(np.nan) - eff_CODs.append(np.nan) - bgs_CODs.append(np.nan) - -rms =[1-i/j for i, j in zip(eff_CODs, inf_CODs)] -outs_CODs = [i+j for i, j in zip(eff_CODs, bgs_CODs)] -balances = [i/j for i, j in zip(outs_CODs, inf_CODs)] - - -# %% -from exposan.metab_mock import * - -sysA, sysB, sysC = create_systems() -dct = globals() -for sys in (sysA, sysB, sysC): dct.update(sys.flowsheet.to_dict()) - -cmps = R1.components -def get_COD(streams): - try: iter(streams) - except: streams = [streams] - return sum((s.mass*cmps.i_COD).sum() for s in streams) - -t = 200 -# t = 20 # shorter t doesn't make sense -# t = 2000 # long t doesn't make a difference (sysB is about to close COD MB, sysC a little off) - -split = 0.9 -R1.split = R2.split = (split, 1-split) - -degassing = 1 -for u in (DM1_c, DM2_c): - for gas in ('H2', 'CH4', 'CO2'): - setattr(u, f'{gas}_degas_efficiency', degassing) -# DM1_c.CO2_degas_efficiency = DM2_c.CO2_degas_efficiency = 0 - - -sysB.simulate(state_reset_hook='reset_cache', t_span=(0,t), method='BDF') -sysC.simulate(state_reset_hook='reset_cache', t_span=(0,t), method='BDF') - -get_bgs = lambda sys: [s for s in sys.products if 'biogas' in s.ID] - -print(f'\nt is {t} d, split is {split}, degassing is {degassing}:') -COD_effB = get_COD(Effluent_B) -COD_effC = get_COD(Effluent_C) -print('sysB eff COD: ', COD_effB) -print('sysC eff COD: ', COD_effC) - -print('sysB biogas COD: ', get_COD(get_bgs(sysB))) -print('sysC biogas COD: ', get_COD(get_bgs(sysC))) - -COD_in = get_COD(sysB.feeds) -assert COD_in == get_COD(sysC.feeds) -print('tot COD in: ', COD_in) -COD_outB = get_COD(sysB.products) -COD_outC = get_COD(sysC.products) -print('sysB tot COD out: ', COD_outB) -print('sysC tot COD out: ', COD_outC) -print('sysB tot COD recovery: ', COD_outB/COD_in) -print('sysC tot COD recovery: ', COD_outC/COD_in) -print('sysB tot COD removal: ', 1-COD_effB/COD_in) -print('sysC tot COD removal: ', 1-COD_effC/COD_in) - -print(f'R1 MB: {get_COD(R1.outs)/get_COD(R1.ins)}') -print(f'R2 MB: {get_COD(R2.outs)/get_COD(R2.ins)}') -print(f'DM1_c MB: {get_COD(DM1_c.outs)/get_COD(DM1_c.ins)}') -print(f'DM2_c MB: {get_COD(DM2_c.outs)/get_COD(DM2_c.ins)}') - - -# %% - -from qsdsan import sanunits as su, WasteStream, processes as pc, System -from exposan.metab_mock import * - -Q = 5 # influent flowrate [m3/d] -T1 = 273.15+35 # temperature [K] -Vl1 = 5 # liquid volume [m^3] -Vg1 = 0.556 # headspace volume [m^3] -split_1 = 0.75 # split ratio to side-stream -tau_1 = 0.021 # degassing membrane retention time [d] - -T2 = 273.15+25 -Vl2 = 75 -Vg2 = 5 -split_2 = 0.75 -tau_2 = 0.021 - -fermenters = ('X_su', 'X_aa', 'X_fa', 'X_c4', 'X_pro') -methanogens = ('X_ac', 'X_h2') -biomass_IDs = (*fermenters, *methanogens) - -adm1 = pc.ADM1() - -inf_b = WasteStream('BreweryWW_B', T=T1) -inf_b.set_flow_by_concentration(Q, concentrations=default_inf_concs, units=('m3/d', 'kg/m3')) -eff_B = WasteStream('Effluent_B', T=T2) -bg1_B = WasteStream('biogas_1B', phase='g') -bg2_B = WasteStream('biogas_2B', phase='g') -AnR1 = su.AnaerobicCSTR('AnR1', ins=inf_b, outs=(bg1_B, ''), - V_liq=Vl1, V_gas=Vg1, T=T1, model=adm1, - retain_cmps=fermenters) -# AnR2 = su.AnaerobicCSTR('AnR2', ins=AnR1-1, outs=(bg2_B, eff_B), -# V_liq=Vl2, V_gas=Vg2, T=T2, model=adm1, -# retain_cmps=methanogens) -AnR1.set_init_conc(**R1_ss_conds) -# AnR2.set_init_conc(**R2_ss_conds) -sysB = System('baseline', path=(AnR1,)) -sysB.set_dynamic_tracker(AnR1, bg1_B) -# sysB = System('baseline', path=(AnR1, AnR2)) -# sysB.set_dynamic_tracker(AnR1, AnR2, bg1_B, bg2_B) - -inf_c = inf_b.copy('BreweryWW_C') -eff_c = WasteStream('Effluent_C', T=T2) -bgm1 = WasteStream('biogas_mem_1', phase='g') -bgm2 = WasteStream('biogas_mem_2', phase='g') -bgh1 = WasteStream('biogas_hsp_1', phase='g') -bgh2 = WasteStream('biogas_hsp_2', phase='g') - -############# sysC unit operation ################# -sc1 = 0.9 -sc2 = 0.9 -R1 = su.AnaerobicCSTR('R1', ins=[inf_c, 'return_1'], - outs=(bgh1, 'sidestream_1', ''), - split=(sc1, 1-sc1), - V_liq=Vl1, V_gas=Vg1, T=T1, model=adm1, - retain_cmps=fermenters) -DM1c = DM1_c = DegassingMembrane('DM1_c', ins=R1-1, outs=(bgm1, 1-R1), tau=tau_1) -# DM1c = DM('DM1_c', ins=R1-1, outs=(bgm1, 1-R1), tau=0.1) - -# R2 = su.AnaerobicCSTR('R2', ins=[R1-2, 'return_2'], -# outs=(bgh2, 'sidestream_2', eff_c), -# split=(sc2, 1-sc2), -# V_liq=Vl2, V_gas=Vg2, T=T2, model=adm1, -# retain_cmps=methanogens) -# DM2c = DM2_c = DegassingMembrane('DM2_c', ins=R2-1, outs=(bgm2, 1-R2), tau=tau_2) -# DM2c = DM('DM2_c', ins=R2-1, outs=(bgm2, 1-R2), tau=0.1) - -R1.set_init_conc(**R1_ss_conds) -# R2.set_init_conc(**R2_ss_conds) - -degassing = 1 -for u in (DM1_c, ): - for gas in ('H2', 'CH4', 'CO2'): - setattr(u, f'{gas}_degas_efficiency', degassing) -# DM1_c.H2_degas_efficiency = DM2_c.H2_degas_efficiency = 1 -# DM1_c.CH4_degas_efficiency = DM2_c.CH4_degas_efficiency = 1 - -sysC = System('combined_METAB', path=(R1, DM1c,), - recycle=(DM1c-1,)) -sysC.set_dynamic_tracker(R1, bgm1, bgh1,) - -# sysC = System('combined_METAB', path=(R1, DM1c, R2, DM2c), -# recycle=(DM1c-1, DM2c-1)) -# sysC.set_dynamic_tracker(R1, R2, bgm1, bgm2, bgh1, bgh2) - - - -cmps = AnR1.components -def get_COD(streams): - try: iter(streams) - except: streams = [streams] - return sum((s.mass*cmps.i_COD).sum() for s in streams) - -t = 200 -# sysB.simulate(state_reset_hook='reset_cache', t_span=(0,t), method='BDF') -sysC.simulate(state_reset_hook='reset_cache', t_span=(0,t), method='BDF') - -#%% -import exposan.metab_mock as mm -sysA, sysB, sysC = mm.create_systems() -# sysC._setup() -# sysC.converge() -cmps = sysC.units[0].components -sysC.set_tolerance(rmol=1e-6) -sysC.simulate(t_span=(0,200), method='BDF', state_reset_hook='reset_cache') -inf_cod = sum([(ws.mass * cmps.i_COD).sum() for ws in sysC.feeds]) -out_cod = sum([(ws.mass * cmps.i_COD).sum() for ws in sysC.products]) -out_cod/inf_cod diff --git a/exposan/metab_mock/graveyard/let23/__init__.py b/exposan/metab_mock/graveyard/let23/__init__.py deleted file mode 100644 index b2016757..00000000 --- a/exposan/metab_mock/graveyard/let23/__init__.py +++ /dev/null @@ -1,34 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Sat Dec 10 13:36:48 2022 - -@author: joy_c -""" - -import os -folder = os.path.dirname(__file__) -data_path = os.path.join(folder, 'data') -results_path = os.path.join(folder, 'results') -figures_path = os.path.join(folder, 'figures') - -# To save simulation results and generated figures -if not os.path.isdir(results_path): os.mkdir(results_path) -if not os.path.isdir(figures_path): os.mkdir(figures_path) -del os - -from . import system -from .system import * - -from . import model -from .model import * - -from . import mock_centralized as mc - -__all__ = ( - 'folder', - 'data_path', - 'results_path', - 'figures_path', - *system.__all__, - *model.__all__, - ) \ No newline at end of file diff --git a/exposan/metab_mock/graveyard/let23/analysis.py b/exposan/metab_mock/graveyard/let23/analysis.py deleted file mode 100644 index 9c4a9c4a..00000000 --- a/exposan/metab_mock/graveyard/let23/analysis.py +++ /dev/null @@ -1,252 +0,0 @@ -# -*- coding: utf-8 -*- -''' -EXPOsan: Exposition of sanitation and resource recovery systems - -This module is developed by: - - Joy Zhang - -This module is under the University of Illinois/NCSA Open Source License. -Please refer to https://github.com/QSD-Group/EXPOsan/blob/main/LICENSE.txt -for license details. -''' - -from qsdsan.utils import load_data, ospath -from exposan.metab_mock.let23 import ( - results_path, - figures_path, - mc, - create_models - ) -import numpy as np, pandas as pd -import matplotlib as mpl, matplotlib.pyplot as plt - -mpl.rcParams['font.sans-serif'] = 'arial' -mpl.rcParams["figure.autolayout"] = False -mpl.rcParams['ytick.minor.visible'] = True - -#%% metric data processing -def load_metrics(seed): - slc = {} - non = {} - for i in 'ABCD': - df = load_data(ospath.join(results_path, f'sys{i}_selective_table_{seed}.xlsx'), - index_col=0, header=[0,1]) - slc[i] = df.iloc[:,23:] - df = load_data(ospath.join(results_path, f'sys{i}_nonselect_table_{seed}.xlsx'), - index_col=0, header=[0,1]) - non[i] = df.iloc[:,23:] - - metrics = {} - iterables = [['selective', 'nonselect'], ['A','B','C','D']] - for col in slc['A'].columns: - l = [slc[i].loc[:,col] for i in 'ABCD'] \ - + [non[i].loc[:,col] for i in 'ABCD'] - metrics[col[1]] = pd.concat(l, axis=1, keys=pd.MultiIndex.from_product(iterables, - names=['group', 'sys'])) - return metrics - -#%% Boxplots of all UA results - -idx = pd.IndexSlice -flierprops = dict(marker='.', markersize=1, markerfacecolor='#90918e', markeredgecolor='#90918e') -meanprops = dict(marker='^', markersize=2.5, markerfacecolor='black', markeredgecolor='black') -medianprops = dict(color='black') - -def plot_dist(df, save_to=None): - fig, axes = plt.subplots(1, 4, sharey=True, figsize=(5, 4)) - for i, sys in enumerate('ABCD'): - ax = axes[i] - data = df.loc[:, idx[:, sys]].to_numpy() - left = data[:,0][~np.isnan(data[:,0])] - right = data[:,1][~np.isnan(data[:,1])] - ax.boxplot([left, right], showmeans=True, labels=['',''], - flierprops=flierprops, meanprops=meanprops, - medianprops=medianprops) - parts = ax.violinplot([left, right], showextrema=False) - for p, c in zip(parts['bodies'], ('#60c1cf', '#F98F60')): - p.set_facecolor(c) - p.set_alpha(1) - ax.tick_params(axis='y', which='both', direction='inout') - ax.set_title(sys) - fig.subplots_adjust(wspace=0) - if save_to: - fig.savefig(ospath.join(figures_path, save_to), dpi=300, - facecolor='white') - else: return fig, axes - -def plot_all_metrics(metrics): - mpl.rcParams['xtick.minor.visible'] = False - for k, df in metrics.items(): - var = k.split(" [")[0] - plot_dist(df, f'{var}.png') - -def plot_E_per_rCOD(metrics): - mpl.rcParams['xtick.minor.visible'] = False - E_use = - metrics['Net operation energy [kW]']/metrics['Hourly COD removal [kgCOD/hr]'] - fig, axes = plot_dist(E_use) - metro_E = mc.E_model() - p5, p25, p50, p75, p95 = np.percentile(metro_E, [5, 25, 50, 75, 95]) - fig, axes = plt.subplots(1, 4, sharey=True, figsize=(5, 4)) - clr = '#9c4b50' - for i, sys in enumerate('ABCD'): - ax = axes[i] - ax.axhspan(p5, p95, color=clr, alpha=0.2, edgecolor=None) - ax.axhspan(p25, p75, color=clr, alpha=0.2) - ax.axhline(y=p50, color=clr, linewidth=0.5) - data = E_use.loc[:, idx[:, sys]].to_numpy() - left = data[:,0][~np.isnan(data[:,0])] - right = data[:,1][~np.isnan(data[:,1])] - ax.boxplot([left, right], showmeans=True, labels=['',''], - flierprops=flierprops, meanprops=meanprops, - medianprops=medianprops) - parts = ax.violinplot([left, right], showextrema=False) - for p, c in zip(parts['bodies'], ('#60c1cf', '#F98F60')): - p.set_facecolor(c) - p.set_alpha(1) - ax.tick_params(axis='y', which='both', direction='inout') - ax.set_ylim(0, 10) - ax.set_title(sys) - fig.subplots_adjust(wspace=0) - fig.savefig(ospath.join(figures_path, 'EprCOD'), dpi=300, - facecolor='white') - - - -#%% SA on sysC (non-selective) net operation energy (MCF) - -import seaborn as sns -from math import ceil -from scipy.stats import kstest - -def plot_cdf_bygroup(x_df, group, save_to): - mpl.rcParams['xtick.minor.visible'] = True - n_param = x_df.shape[1] - ncol = 3 - nrow = ceil(n_param/ncol) - fig, axes = plt.subplots(nrow, ncol, sharey=True, figsize=(ncol*2, nrow*2.5)) - for col, ax in zip(x_df, axes.ravel()): - sns.ecdfplot(data=x_df, x=col, hue=group, - palette=['#a280b9', '#f3c354'], - legend=False, ax=ax) - ax.tick_params(axis='both', which='both', direction='inout') - ax.set_xlabel('') - ax.set_ylabel('') - fig.subplots_adjust(hspace=0.4, wspace=0.05, bottom=0.2) - fig.savefig(ospath.join(figures_path, save_to), dpi=300, facecolor='white') - - -def KStest(seed, yname, threshold, sys='sysC', scenario='nonselect', plot_save_to='cdf.png'): - df = load_data(ospath.join(results_path, f'{sys}_{scenario}_table_{seed}.xlsx'), - index_col=0, header=[0,1]) - ys = df.xs(key=yname, axis=1, level='Feature').to_numpy() - highx = df.loc[ys > threshold, idx[['Anaerobic CSTR-R1C', 'ADM1'], :]] - lowx = df.loc[ys <= threshold, idx[['Anaerobic CSTR-R1C', 'ADM1'], :]] - stats = {} - sig_cols = [] - for col in highx.columns: - x = col[1].split(' [')[0] - D, p = kstest(highx.loc[:, col], lowx.loc[:, col]) - stats[x] = (D, p) - if p < 0.05: sig_cols.append(col) - group = ['high']*highx.shape[0] + ['low']*lowx.shape[0] - xs = pd.concat([highx, lowx])[sig_cols] - plot_cdf_bygroup(xs, group, plot_save_to) - return stats, sig_cols - -#%% Baseline energy breakdown -path = ospath.join(results_path, 'E_breakdown.xlsx') - -def baseline_E_breakdown(): - smdls = create_models(selective=True) - nmdls = create_models(selective=False) - - outs = [] - for mdl in smdls+nmdls: - mdl.system.simulate(t_span=(0,400), method='BDF', state_reset_hook='reset_cache') - outs.append([m.getter() for m in mdl.metrics]) - cols = [f'{m.name} [{m.units}]' for m in mdl.metrics] - iterables = [['selective', 'nonselect'], ['A','B','C','D']] - rows = pd.MultiIndex.from_product(iterables, names=['group', 'sys']) - outs = pd.DataFrame(outs, columns=cols, index=rows) - outs.to_excel(path) - -df = load_data(path, sheet=1, index_col=[0,1]) -viridis = plt.cm.get_cmap('viridis', 5).colors -def stacked_bar(df, save_to=None): - mpl.rcParams['xtick.minor.visible'] = False - fig, axes = plt.subplots(1, 4, sharey=True, figsize=(5, 4)) - x = [0.3, 0.7] - for i, sys in enumerate('ABCD'): - ax = axes[i] - ax.axhline(y=0, color='black', linewidth=0.5) - data = df.loc[idx[:, sys],:].to_numpy().T - yp = np.zeros(2) - yn = yp.copy() - for j, y in enumerate(data): - y_offset = (y>=0)*yp + (y<0)*yn - ax.bar(x, y, width=0.2, bottom=y_offset, color=viridis[j], - tick_label='') - yp += (y>=0) * y - yn += (y<0) * y - ax.tick_params(axis='y', which='both', direction='inout') - ax.set_xlim(0,1) - ax.set_title(sys) - fig.subplots_adjust(wspace=0) - if save_to: - fig.savefig(ospath.join(figures_path, save_to), dpi=300, - transparent=True) - else: return fig, axes - -#%% -def compare_metrics(metrics): - compare = {} - iterables = [['selective', 'nonselect'], ['B-A','C-A','D-C']] - keys = pd.MultiIndex.from_product(iterables, names=['group', 'pair']) - for k, df in metrics.items(): - l = [] - for g in iterables[0]: - l.append(df[g].B - df[g].A) - l.append(df[g].C - df[g].A) - l.append(df[g].D - df[g].C) - # l.append((df[g].B/df[g].A-1)*100) - # l.append((df[g].C/df[g].A-1)*100) - # l.append((df[g].D/df[g].C-1)*100) - compare[k] = pd.concat(l, axis=1, keys=keys) - return compare - -def plot_compare(df, save_to=None): - fig, axes = plt.subplots(1, 2, sharey=True, figsize=(3, 4)) - for g, i, c in zip(['selective', 'nonselect'],(0,1), ('#60c1cf', '#F98F60')): - ax = axes[i] - ax.axhline(y=0, color='black', linewidth=0.5) - cols = [col[~np.isnan(col)] for col in df[g].to_numpy().T] - ax.boxplot(cols, showmeans=True, labels=['B-A','C-A','D-C'], - flierprops=flierprops, meanprops=meanprops, - medianprops=medianprops) - parts = ax.violinplot(cols, showextrema=False) - for p in parts['bodies']: - p.set_facecolor(c) - p.set_alpha(1) - ax.tick_params(axis='y', which='both', direction='inout') - fig.subplots_adjust(wspace=0) - if save_to: - fig.savefig(ospath.join(figures_path, f'diff/{save_to}'), dpi=300, - facecolor='white', bbox_inches='tight') - else: return fig, axes - -def plot_all_compare(compare): - mpl.rcParams['xtick.minor.visible'] = False - for k, df in compare.items(): - var = k.split(" [")[0] - plot_compare(df, f'{var}.png') - - -#%% -if __name__ == '__main__': - seed = 123 - metrics = load_metrics(seed) - # plot_all_metrics(metrics) - # plot_E_per_rCOD(metrics) - # yname = 'Net operation energy [kW]' - # stats, sig_xs = KStest(seed, yname, threshold=-1.15) diff --git a/exposan/metab_mock/graveyard/let23/mock_centralized.py b/exposan/metab_mock/graveyard/let23/mock_centralized.py deleted file mode 100644 index e4b9da66..00000000 --- a/exposan/metab_mock/graveyard/let23/mock_centralized.py +++ /dev/null @@ -1,45 +0,0 @@ -# -*- coding: utf-8 -*- -''' -EXPOsan: Exposition of sanitation and resource recovery systems - -This module is developed by: - - Joy Zhang - -This module is under the University of Illinois/NCSA Open Source License. -Please refer to https://github.com/QSD-Group/EXPOsan/blob/main/LICENSE.txt -for license details. -''' -from qsdsan.utils import auom -import chaospy.distributions as shape, numpy as np -import chaospy - -N_blower = 7 -lhp = 3000 -hhp = 5500 - - -def E_per_rCOD(n_lhp, Q, BOD_in, rBOD_pri, BOD2COD_pre2nd, BOD_eff, BOD2COD_eff): - avg_power = (n_lhp * lhp + (N_blower - n_lhp) * hhp) \ - * auom('hp').conversion_factor('kW') - Q = Q * auom('mgd').conversion_factor('L/hr') - pre2_BOD = BOD_in * (1-rBOD_pri) - pre2_COD = pre2_BOD / BOD2COD_pre2nd - eff_COD = BOD_eff / BOD2COD_eff - rCOD = (pre2_COD - eff_COD) * auom('mg/L').conversion_factor('kg/L') - return avg_power / (rCOD * Q) # kWh/kg COD removal - - -def E_model(): - RVs = [] - RVs.append(shape.Binomial(N_blower, 0.5)) # distribution of the number of low-hp blowers - RVs.append(shape.Triangle(170, 185, 250)) # distribution of influent Q [MGD] - RVs.append(shape.Triangle(110, 190, 350)) # untreaded municipal WW BOD [mg/L] - RVs.append(shape.Uniform(0.5, 0.6)) # primary treatment BOD removal rate - RVs.append(shape.Uniform(0.4, 0.6)) # BOD/COD after primary treatment - RVs.append(shape.Triangle(0, 5, 5)) # Effluent BOD - RVs.append(shape.Uniform(0.1, 0.3)) # BOD/COD after secondary treatment - params = chaospy.J(*RVs) - samples = params.sample(size=400, rule='latin_hypercube').T - Es = np.array([E_per_rCOD(*row) for row in samples]) - return Es \ No newline at end of file diff --git a/exposan/metab_mock/graveyard/let23/model.py b/exposan/metab_mock/graveyard/let23/model.py deleted file mode 100644 index 64d14ee2..00000000 --- a/exposan/metab_mock/graveyard/let23/model.py +++ /dev/null @@ -1,331 +0,0 @@ -# -*- coding: utf-8 -*- -''' -EXPOsan: Exposition of sanitation and resource recovery systems - -This module is developed by: - - Joy Zhang - -This module is under the University of Illinois/NCSA Open Source License. -Please refer to https://github.com/QSD-Group/EXPOsan/blob/main/LICENSE.txt -for license details. -''' - -import qsdsan as qs, qsdsan.sanunits as su, numpy as np, os -# from warnings import warn -from chaospy import distributions as shape -from qsdsan.utils import ospath, time_printer, MethodSetter, auom -from exposan.metab_mock.let23 import create_systems, results_path -# import pandas as pd - -__all__ = ('create_models', - 'run_model', - 'run_uncertainty') -#%% -ks = { - 'disintegration': (0.5, 3), # baseline, variation factor - 'hydrolysis_carbs': (10, 1), - 'hydrolysis_proteins': (10, 1), - 'hydrolysis_lipids': (10, 3), - 'uptake_LCFA': (6, 3), - 'uptake_propionate': (13, 1), - 'uptake_acetate': (8, 1), - } -# bs = (0.02, 1) -# fb = (1, 1, 2) # baseline, lower bound, upper bound - -Ks = { - 'uptake_LCFA': (0.4, 3), - 'uptake_c4': (0.2, 3), - 'uptake_propionate': (0.1, 1), - 'uptake_acetate': (0.15, 1), - 'uptake_h2': (7e-6, 1) - } - -# KIs_h2 = { -# 'uptake_propionate': (3.5e-6, 0.3), -# } - -# KI_nh3 = (0.0018, 0.3) -# pH_ULs = { -# 'uptake_acetate': (7, 0.3), -# 'uptake_h2': (6, 0.5) -# } -# pH_LLs = { -# 'uptake_acetate': (6, 0.3), -# 'uptake_h2': (5, 0.3) -# } - -def get_uniform_from_fvar(b, f): - u = b * (1+f) - if f >= 1: l = b*0.01 - else: l = b * (1-f) - return shape.Uniform(l, u) - -def add_params(model): - param = model.parameter - R1, R2 = units = [u for u in model.system.units if isinstance(u, su.AnaerobicCSTR)] - - #****** technological parameters ******* - b = 0.1 - D = shape.Uniform(0.1, 0.5) - @param(name='Headspace_P', element=R1, kind='coupled', units='atm', - baseline=b, distribution=D) - def set_hsp_P(P): - if R1.fixed_headspace_P: R1.headspace_P = P*1.01325 - - b = 400 - D = shape.Uniform(1, 400) - @param(name='Recirculation_rate', element=R1, kind='coupled', units='', - baseline=b, distribution=D) - def set_recir(r): - if R1.split is not None: R1.split = [r, 1] - - b = 0.95 - D = shape.Uniform(0.9, 1.0) - @param(name='Biomass_retention', element=R1, kind='coupled', units='', - baseline=b, distribution=D) - def set_f_retain(f): - for u in units: - u._f_retain = (u._f_retain > 0) * f - - #****** ADM1 parameters ********* - adm1 = R1.model - for p, v in ks.items(): - b, f = v - D = get_uniform_from_fvar(b, f) - setter = MethodSetter(adm1, 'set_rate_constant', process=p) - name = f'mu_{p}' - param(setter=setter, name=name, element=adm1, kind='coupled', units='d^(-1)', - baseline=b, distribution=D) - - b = 0.02 - D = get_uniform_from_fvar(b, 1) - @param(name='Slow_decay_b', element=adm1, kind='coupled', - units='d^(-1)', baseline=b, distribution=D) - def set_slow_b(k): - adm1.rate_function._params['rate_constants'][12:19] = k - - b = 1 - D = shape.Uniform(1, 2) - idx = adm1.indices(['decay_Xsu', 'decay_Xaa', 'decay_Xac']) - @param(name='Fast_to_slow_b_ratio', element=adm1, kind='coupled', units='', - baseline=b, distribution=D) - def set_b_ratio(r): - adm1.rate_function._params['rate_constants'][idx] *= r - - for p, v in Ks.items(): - b, f = v - D = get_uniform_from_fvar(b, f) - if p.endswith('_c4'): - @param(name='K_uptake_c4', element=adm1, kind='coupled', units='kgCOD/m3', - baseline=b, distribution=D) - def set_Ks_c4(K): - adm1.set_half_sat_K(K, process='uptake_butyrate') - adm1.set_half_sat_K(K, process='uptake_valerate') - else: - setter = MethodSetter(adm1, 'set_half_sat_K', process=p) - param(setter=setter, name=f'K_{p}', element=adm1, kind='coupled', - units='kgCOD/m3', baseline=b, distribution=D) - - b = 3.5e-6 - D = get_uniform_from_fvar(b, 0.3) - @param(name='KI_h2_uptake_propionate', element=adm1, kind='coupled', - units='kgCOD/m3', baseline=b, distribution=D) - def set_KIh2_pro(KI): - adm1.set_h2_inhibit_K(KI, process='uptake_propionate') - - b = 0.0018 - D = get_uniform_from_fvar(b, 0.3) - @param(name='KI_nh3', element=adm1, kind='coupled', units='M', - baseline=b, distribution=D) - def set_KI_nh3(KI): - adm1.set_KI_nh3(KI) - - b = 7 - D = get_uniform_from_fvar(b, 0.3) - @param(name='pH_UL_uptake_acetate',element=adm1, kind='coupled', units='', - baseline=b, distribution=D) - def set_pH_UL_ac(ul): - adm1.rate_function._params['pH_ULs'][-2] = ul - - b = 6 - D = get_uniform_from_fvar(b, 0.5) - @param(name='pH_UL_uptake_h2',element=adm1, kind='coupled', units='', - baseline=b, distribution=D) - def set_pH_UL_h2(ul): - adm1.rate_function._params['pH_ULs'][-1] = ul - - b = 6 - D = get_uniform_from_fvar(b, 0.3) - @param(name='pH_LL_uptake_acetate',element=adm1, kind='coupled', units='', - baseline=b, distribution=D) - def set_pH_LL_ac(ll): - ULs = adm1.rate_function._params['pH_ULs'] - LLs = adm1.rate_function._params['pH_LLs'] - ul = ULs[-2] - if ll < ul: LLs[-2] = ll - elif ll == ul: - ULs[-2] = ul + 0.5 - LLs[-2] = ll - 0.5 - else: - ULs[-2] = ll - LLs[-2] = ul - - b = 5 - D = get_uniform_from_fvar(b, 0.3) - @param(name='pH_LL_uptake_h2',element=adm1, kind='coupled', units='', - baseline=b, distribution=D) - def set_pH_LL_h2(ll): - ULs = adm1.rate_function._params['pH_ULs'] - LLs = adm1.rate_function._params['pH_LLs'] - ul = ULs[-1] - if ll < ul: LLs[-1] = ll - elif ll == ul: - ULs[-1] = ul + 0.5 - LLs[-1] = ll - 0.5 - else: - ULs[-1] = ll - LLs[-1] = ul - -def add_metrics(model): - metric = model.metric - sys = model.system - inf, = sys.feeds - eff, = [ws for ws in sys.products if ws.ID.startswith('Effluent')] - bgs = [ws for ws in sys.products if ws.ID.startswith('Biogas')] - pumps = [u for u in sys.units if isinstance(u, su.Pump)] - reactors = [u for u in sys.units if isinstance(u, su.AnaerobicCSTR)] - - cmps = eff.components - S_h2_i_mass = cmps.S_h2.i_mass - S_ch4_i_mass = cmps.S_ch4.i_mass - cmps_i_COD = cmps.i_COD - exclude_gas = cmps.s + cmps.x - substrate_IDs = ('S_su', 'S_aa', 'S_fa', 'S_va', 'S_bu', 'S_pro', 'S_ac', - 'X_c', 'X_ch', 'X_pr', 'X_li') - - @metric(name='Soluble COD', units='mg/L', element='Effluent') - def get_sCOD(): - return eff.composite('COD', particle_size='s') - - @metric(name='Particulate COD', units='mg/L', element='Effluent') - def get_xCOD(): - return eff.composite('COD', particle_size='x') - - @metric(name='Substrate COD', units='mg/L', element='Effluent') - def get_subCOD(): - return eff.composite('COD', subgroup=substrate_IDs) - - @metric(name='Total COD removal', units='%', element='System') - def get_rCOD(): - return (1 - sum(eff.mass*cmps_i_COD*exclude_gas)/sum(inf.mass*cmps_i_COD*exclude_gas))*100 - - @metric(name='Hourly COD removal', units='kgCOD/hr', element='System') - def get_rCOD_hr(): - return sum(inf.mass*cmps_i_COD*exclude_gas) - sum(eff.mass*cmps_i_COD*exclude_gas) # kg/hr removal - - @metric(name='H2 production', units='kg/d', element='Biogas') - def get_QH2(): - return sum([bg.imass['S_h2'] for bg in bgs])*S_h2_i_mass*24 - - @metric(name='CH4 production', units='kg/d', element='Biogas') - def get_QCH4(): - return sum([bg.imass['S_ch4'] for bg in bgs])*S_ch4_i_mass*24 - - @metric(name='H2 yield', units='kg/kgCOD-removed', element='Biogas') - def get_yH2(): - rCOD = sum(inf.mass*cmps_i_COD*exclude_gas) - sum(eff.mass*cmps_i_COD*exclude_gas) # kg/hr removal - return sum([bg.imass['S_h2'] for bg in bgs])*S_h2_i_mass/rCOD - - @metric(name='CH4 yield', units='kg/kgCOD-removed', element='Biogas') - def get_yCH4(): - rCOD = sum(inf.mass*cmps_i_COD*exclude_gas) - sum(eff.mass*cmps_i_COD*exclude_gas) # kg/hr removal - return sum([bg.imass['S_ch4'] for bg in bgs])*S_ch4_i_mass/rCOD - - @metric(name='Pumping electricity', units='kW', element='Energy') - def get_Epump(): - return sum(p.parallel['self']*p.power_utility.consumption for p in pumps) - - @metric(name='Heat', units='MJ/hr', element='Energy') - def get_Eheat(): - return reactors[0].heat_utilities[0].duty/1e3 # kJ/hr to MJ/hr - - @metric(name='Biogas energy', units='MJ/hr', element='Energy') - def get_E_biogas(): - # kmol/hr * J/mol = kJ/hr - return sum(sum(bg.mass*cmps.i_mass/cmps.chem_MW*cmps.LHV) for bg in bgs)/1e3 - - @metric(name='H2 energy', units='MJ/hr', element='Energy') - def get_E_h2(): - return sum(bg.imass['S_h2'] for bg in bgs)*S_h2_i_mass/cmps.S_h2.chem_MW*cmps.S_h2.LHV/1e3 - - @metric(name='CH4 energy', units='MJ/hr', element='Energy') - def get_E_ch4(): - return sum(bg.imass['S_ch4'] for bg in bgs)*S_ch4_i_mass/cmps.S_ch4.chem_MW*cmps.S_ch4.LHV/1e3 - - @metric(name='H2O energy', units='MJ/hr', element='Energy') - def get_E_h2o(): - return sum(bg.imass['H2O'] for bg in bgs)/cmps.H2O.MW*cmps.H2O.LHV/1e3 - - conv = auom('kJ/hr').conversion_factor('kW') - @metric(name='Net operation energy', units='kW', element='Energy') - def get_net_E(): - consume = sum(p.parallel['self']*p.power_utility.consumption for p in pumps) \ - + reactors[0].heat_utilities[0].duty * conv - produce = sum(sum(bg.mass*cmps.i_mass/cmps.chem_MW*cmps.LHV) \ - for bg in bgs) * conv * 0.6 # assume 60% CHP efficiency - return produce - consume - -#%% -def create_models(systems=None, **kwargs): - systems = systems or create_systems(**kwargs) - models = [] - for sys in systems: - mdl = qs.Model(sys, exception_hook='warn') - add_params(mdl) - add_metrics(mdl) - models.append(mdl) - return models - -@time_printer -def run_model(model, N, T, t_step, method='BDF', - metrics_path='', timeseries_path='', prefix='', - rule='L', seed=None, pickle=False): - if seed: np.random.seed(seed) - sys_ID = model._system.ID - samples = model.sample(N=N, rule=rule) - model.load_samples(samples) - t_span = (0, T) - t_eval = np.arange(0, T+t_step, t_step) - mpath = metrics_path or ospath.join(results_path, f'{sys_ID}_{prefix}_table_{seed}.xlsx') - if timeseries_path: tpath = timeseries_path - else: - folder = ospath.join(results_path, f'{sys_ID}_{prefix}_ytdata_{seed}') - os.mkdir(folder) - tpath = ospath.join(folder, 'state.npy') - model.evaluate( - state_reset_hook='reset_cache', - t_span=t_span, - t_eval=t_eval, - method=method, - export_state_to=tpath - ) - model.table.to_excel(mpath) - -def run_uncertainty(seed, N, models=None, T=400, t_step=5, method='BDF', selective=True): - mdls = models or create_models(selective=selective) - if selective: prefix = 'selective' - else: prefix = 'nonselect' - for mdl in mdls: - msg = f'Model {mdl.system.ID[-1]} ({prefix})' - print(f'\n{msg}\n{"-"*len(msg)}') - print(f'Seed = {seed} \nN = {N} \nTime span 0-{T}d \n') - run_model(mdl, N=N, T=T, t_step=t_step, prefix=prefix, seed=seed) - # return mdls - -#%% -if __name__ == '__main__': - # mdls = create_models(which='D') - # run_uncertainty(123, 10, models=mdls) - run_uncertainty(123, 1000, selective=False) \ No newline at end of file diff --git a/exposan/metab_mock/graveyard/let23/system.py b/exposan/metab_mock/graveyard/let23/system.py deleted file mode 100644 index 658352d9..00000000 --- a/exposan/metab_mock/graveyard/let23/system.py +++ /dev/null @@ -1,300 +0,0 @@ -# -*- coding: utf-8 -*- -''' -EXPOsan: Exposition of sanitation and resource recovery systems - -This module is developed by: - - Joy Zhang - -This module is under the University of Illinois/NCSA Open Source License. -Please refer to https://github.com/QSD-Group/EXPOsan/blob/main/LICENSE.txt -for license details. -''' - -import numpy as np, qsdsan as qs -from qsdsan import sanunits as su, processes as pc, WasteStream, System -from qsdsan.utils import ExogenousDynamicVariable as EDV, auom -# from chemicals.elements import molecular_weight as get_mw -from exposan.metab_mock import DegassingMembrane as DM, rhos_adm1_ph_ctrl - -__all__ = ('create_systems',) - -#%% Global Variables - -from exposan.metab_mock.systems import ( - Q, T1, Vl1, Vg1, T2, Vl2, Vg2, - fermenters, methanogens, biomass_IDs, #vfa_IDs, - default_inf_concs, R1_ss_conds, R2_ss_conds, - # yields_bl, mus_bl, Ks_bl, - ) - -# reactor water depths [m] - -h1 = 4.5 -h2 = 6 - -q = Q*auom('m3/d').conversion_factor('ft3/s') -L = 5 # pipe length [ft] -D = 3*auom('inch').conversion_factor('ft') -C = 110 # the Hazen-Williams coefficient -hf = 3.02*L*(4*q/C/np.pi)*D**(-4.87)*auom('ft').conversion_factor('m') # friction head -dP1 = 1000 * 9.8 * (h1+hf*2) -dP2 = 1000 * 9.8 * (h2-h1+hf*2) -dPs = 1000 * 9.8 * hf*2 -dPm = 6*auom('psi').conversion_factor('Pa') # Hofi membrane water pressure drop at median flowrate through lumen - -ph1 = 5.8 -ph2 = 7.2 - -Temp1 = EDV('T1', function=lambda t: T1) -Temp2 = EDV('T2', function=lambda t: T2) -pH1 = EDV('pH1', function=lambda t: ph1) -pH2 = EDV('pH2', function=lambda t: ph2) - -#%% Systems - -def create_systems(flowsheet_A=None, flowsheet_B=None, flowsheet_C=None, flowsheet_D=None, - inf_concs={}, R1_init_conds={}, R2_init_conds={}, which=None, - selective=True): - - which = which or ('A', 'B', 'C', 'D') - if isinstance(which, str): which = (which,) - - if selective: - R1_retain = fermenters - R2_retain = methanogens - else: R1_retain = R2_retain = biomass_IDs - pc.create_adm1_cmps() - adm1 = pc.ADM1() - dyn_params = adm1.rate_function.params.copy() - adm1.set_rate_function(rhos_adm1_ph_ctrl) - adm1.rate_function._params = dyn_params - inf_concs = inf_concs or default_inf_concs.copy() - brewery_ww = WasteStream('BreweryWW_A') - brewery_ww.set_flow_by_concentration(Q, concentrations=inf_concs, units=('m3/d', 'kg/m3')) - R1_init_conds = R1_init_conds or R1_ss_conds - R2_init_conds = R2_init_conds or R2_ss_conds - systems = [] - - if 'A' in which: - flowsheet_A = flowsheet_A or qs.Flowsheet('A') - qs.main_flowsheet.set_flowsheet(flowsheet_A) - flowsheet_A.stream.register(brewery_ww.ID, brewery_ww) - - effA = WasteStream('Effluent_A', T=T2) - bg1A = WasteStream('Biogas_1A', phase='g') - bg2A = WasteStream('Biogas_2A', phase='g') - - # adm1A = pc.ADM1() - # dyn_params = adm1A.rate_function.params.copy() - # adm1A.set_rate_function(rhos_adm1_ph_ctrl) - # adm1A.rate_function._params = dyn_params - P1A = su.Pump('P1A', ins=brewery_ww, dP_design=dP1, - init_with='WasteStream', isdynamic=True) - R1A = su.AnaerobicCSTR('R1A', ins=P1A-0, outs=(bg1A, ''), - V_liq=Vl1, V_gas=Vg1, T=T1, model=adm1, - retain_cmps=R1_retain, - exogenous_vars=(Temp1, pH1)) - P2A = su.Pump('P2A', ins=R1A-1, dP_design=dP2, - init_with='WasteStream', isdynamic=True) - R2A = su.AnaerobicCSTR('R2A', ins=P2A-0, outs=(bg2A, effA), - V_liq=Vl2, V_gas=Vg2, T=T2, model=adm1, - retain_cmps=R2_retain, - exogenous_vars=(Temp2, pH2)) - - R1A.set_init_conc(**R1_init_conds) - R2A.set_init_conc(**R2_init_conds) - - sysA = System('sysA', path=(P1A, R1A, P2A, R2A)) - sysA.set_dynamic_tracker(R1A, R2A, bg1A, bg2A, effA) - - systems.append(sysA) - - if 'B' in which: - flowsheet_B = flowsheet_B or qs.Flowsheet('B') - qs.main_flowsheet.set_flowsheet(flowsheet_B) - - infB = brewery_ww.copy('BreweryWW_B') - effB = WasteStream('Effluent_B', T=T2) - bg1B = WasteStream('Biogas_1B', phase='g') - bg2B = WasteStream('Biogas_2B', phase='g') - - # adm1B = pc.ADM1() - # dyn_params = adm1B.rate_function.params.copy() - # adm1B.set_rate_function(rhos_adm1_ph_ctrl) - # adm1B.rate_function._params = dyn_params - P1B = su.Pump('P1B', ins=infB, dP_design=dP1, - init_with='WasteStream', isdynamic=True) - R1B = su.AnaerobicCSTR('R1B', ins=P1B-0, outs=('', ''), - V_liq=Vl1, V_gas=Vg1, T=T1, model=adm1, - fixed_headspace_P=True, - headspace_P=0.1013, - retain_cmps=R1_retain, - exogenous_vars=(Temp1, pH1)) - VP1B = su.Pump('VP1B', ins=R1B-0, outs=bg1B, P=101325) - P2B = su.Pump('P2B', ins=R1B-1, dP_design=dP2, - init_with='WasteStream', isdynamic=True) - R2B = su.AnaerobicCSTR('R2B', ins=P2B-0, outs=(bg2B, effB), - V_liq=Vl2, V_gas=Vg2, T=T2, model=adm1, - retain_cmps=R2_retain, - exogenous_vars=(Temp2, pH2)) - - R1B.set_init_conc(**R1_init_conds) - R2B.set_init_conc(**R2_init_conds) - - sysB = System('sysB', path=(P1B, R1B, VP1B, P2B, R2B)) - sysB.set_dynamic_tracker(R1B, R2B, bg1B, bg2B, effB) - - systems.append(sysB) - - if 'C' in which: - flowsheet_C = flowsheet_C or qs.Flowsheet('C') - qs.main_flowsheet.set_flowsheet(flowsheet_C) - - infC = brewery_ww.copy('BreweryWW_C') - effC = WasteStream('Effluent_C', T=T2) - bg1C = WasteStream('Biogas_1C', phase='g') - bgh2C = WasteStream('Biogas_hsp_2C', phase='g') - bgm2C = WasteStream('Biogas_mem_2C', phase='g') - - # adm1C = pc.ADM1() - # dyn_params = adm1C.rate_function.params.copy() - # adm1C.set_rate_function(rhos_adm1_ph_ctrl) - # adm1C.rate_function._params = dyn_params - P1C = su.Pump('P1C', ins=infC, dP_design=dP1, - init_with='WasteStream', isdynamic=True) - R1C = su.AnaerobicCSTR('R1C', ins=P1C-0, outs=(bg1C, ''), - V_liq=Vl1, V_gas=Vg1, T=T1, model=adm1, - retain_cmps=R1_retain, - exogenous_vars=(Temp1, pH1)) - P2C = su.Pump('P2C', ins=R1C-1, dP_design=dP2, - init_with='WasteStream', isdynamic=True) - R2C = su.AnaerobicCSTR('R2C', ins=P2C-0, outs=(bgh2C, ''), - V_liq=Vl2, V_gas=Vg2, T=T2, model=adm1, - retain_cmps=R2_retain, - exogenous_vars=(Temp2, pH2)) - P2sC = su.Pump('P2sC', ins=R2C-1, dP_design=dPm+dPs, - init_with='WasteStream', isdynamic=True) - DM2C = DM('DM2C', ins=P2sC-0, outs=('', effC)) - VP2C = su.Pump('VP2C', ins=DM2C-0, outs=bgm2C, dP_design=101325) - - R1C.set_init_conc(**R1_init_conds) - R2C.set_init_conc(**R2_init_conds) - - sysC = System('sysC', path=(P1C, R1C, P2C, R2C, P2sC, DM2C, VP2C)) - sysC.set_dynamic_tracker(R1C, R2C, bg1C, bgh2C, bgm2C, effC) - - systems.append(sysC) - - if 'D' in which: - flowsheet_D = flowsheet_D or qs.Flowsheet('D') - qs.main_flowsheet.set_flowsheet(flowsheet_D) - - infD = brewery_ww.copy('BreweryWW_D') - effD = WasteStream('Effluent_D', T=T2) - bgh1D = WasteStream('Biogas_hsp_1D', phase='g') - bgm1D = WasteStream('Biogas_mem_1D', phase='g') - bgh2D = WasteStream('Biogas_hsp_2D', phase='g') - bgm2D = WasteStream('Biogas_mem_2D', phase='g') - - # adm1D = pc.ADM1() - # dyn_params = adm1D.rate_function.params.copy() - # adm1D.set_rate_function(rhos_adm1_ph_ctrl) - # adm1D.rate_function._params = dyn_params - P1D = su.Pump('P1D', ins=infD, dP_design=dP1, - init_with='WasteStream', isdynamic=True) - R1D = su.AnaerobicCSTR('R1D', ins=[P1D-0, 'Return_1D'], - outs=(bgh1D, 'Sidestream_1D', ''), - split=[400, 1], - V_liq=Vl1, V_gas=Vg1, T=T1, model=adm1, - retain_cmps=R1_retain, - exogenous_vars=(Temp1, pH1)) - P1sD = su.Pump('P1sD', ins=R1D-1, dP_design=dPm+dPs, - init_with='WasteStream', isdynamic=True) - DM1D = DM('DM1D', ins=P1sD-0, outs=('', 1-R1D)) - VP1D = su.Pump('VP1D', ins=DM1D-0, outs=bgm1D, dP_design=101325) - P2D = su.Pump('P2D', ins=R1D-2, dP_design=dP2, - init_with='WasteStream', isdynamic=True) - R2D = su.AnaerobicCSTR('R2D', ins=P2D-0, outs=(bgh2D, ''), - V_liq=Vl2, V_gas=Vg2, T=T2, model=adm1, - retain_cmps=R2_retain, - exogenous_vars=(Temp2, pH2)) - P2sD = su.Pump('P2sD', ins=R2D-1, dP_design=dPm+dPs, - init_with='WasteStream', isdynamic=True) - DM2D = DM('DM2D', ins=P2sD-0, outs=('', effD)) - VP2D = su.Pump('VP2D', ins=DM2D-0, outs=bgm2D, dP_design=101325) - - R1D.set_init_conc(**R1_init_conds) - R2D.set_init_conc(**R2_init_conds) - - sysD = System('sysD', path=(P1D, R1D, P1sD, DM1D, VP1D, - P2D, R2D, P2sD, DM2D, VP2D), - recycle=(DM1D-1,)) - sysD.set_dynamic_tracker(R1D, R2D, bgh1D, bgm1D, bgh2D, bgm2D, effD) - - systems.append(sysD) - - return systems - - -#%% -if __name__ == '__main__': - sysA, sysB, sysC, sysD = systems = create_systems() - for sys in systems: - sys.simulate( - t_span=(0,400), - state_reset_hook='reset_cache', - method='BDF' - ) - sys.diagram() - au = sysA.flowsheet.unit - bu = sysB.flowsheet.unit - cu = sysC.flowsheet.unit - du = sysD.flowsheet.unit - -#%% -# cmps = qs.get_thermo().chemicals -# #!!! should we include dissolved gas in COD removal calculation? -# def tCOD(streams): -# # return sum((s.mass*cmps.i_COD).sum() for s in streams) -# return sum(s.composite('COD', flow=True) for s in streams) - -# def lcod(streams): -# # return sum((s.mass*cmps.i_COD).sum() for s in streams if s.phase == 'l') -# return sum(s.composite('COD', flow=True) for s in streams if s.phase == 'l') - -# inf_cod = [tCOD(sys.feeds) for sys in systems] -# out_cod = [tCOD(sys.products) for sys in systems] -# balance = [out/inf for inf, out in zip(inf_cod, out_cod)] - -# eff_cod = [lcod(sys.products) for sys in systems] -# rcod = [1-eff/inf for inf, eff in zip(inf_cod, eff_cod)] - -# rcod = [1-eff/inf for inf, eff in zip(out_cod, eff_cod)] - -#%% -# adm1 = du.R1D.model -# sto_h2 = adm1.stoichiometry.S_h2.to_numpy() -# sto_ch4 = adm1.stoichiometry.S_ch4.to_numpy() -# sto_ac = adm1.stoichiometry.S_ac.to_numpy() - -# def r(unit, sto): -# react = sum(unit._tempstate['rhos'] * sto[4:12]) -# transfer = sum(unit._tempstate['gas_transfer'] * sto[-3:]) -# return react, transfer, react+transfer - -# r1h2 = [] -# r1ac = [] -# r1ch4 = [] - -# r2h2 = [] -# r2ac = [] -# r2ch4 = [] -# for sys in systems: -# R1, R2 = [u for u in sys.units if isinstance(u, su.AnaerobicCSTR)] -# r1h2.append(r(R1, sto_h2)) -# r1ac.append(r(R1, sto_ac)) -# r1ch4.append(r(R1, sto_ch4)) -# r2h2.append(r(R2, sto_h2)) -# r2ac.append(r(R2, sto_ac)) -# r2ch4.append(r(R2, sto_ch4)) diff --git a/exposan/metab_mock/graveyard/models.py b/exposan/metab_mock/graveyard/models.py deleted file mode 100644 index 2b9ec082..00000000 --- a/exposan/metab_mock/graveyard/models.py +++ /dev/null @@ -1,367 +0,0 @@ -# -*- coding: utf-8 -*- -''' -EXPOsan: Exposition of sanitation and resource recovery systems - -This module is developed by: - - Joy Zhang - -This module is under the University of Illinois/NCSA Open Source License. -Please refer to https://github.com/QSD-Group/EXPOsan/blob/main/LICENSE.txt -for license details. -''' - -import qsdsan as qs -# from warnings import warn -from chaospy import distributions as shape -from qsdsan.utils import ospath, time_printer, get_SRT, AttrFuncSetter -from exposan.metab_mock import systems as s, results_path, biomass_IDs, vfa_IDs -from biosteam.evaluation._utils import var_indices -# import pandas as pd -import numpy as np -import os - -__all__ = ('create_modelA', - 'create_modelB', - 'create_modelC', - 'create_modelD', - 'create_ss_model', - 'run_model', - 'run_modelB', - 'run_ss_model' - ) -#%% -Ys_bl, mus_bl, Ks_bl = s.yields_bl, s.mus_bl, s.Ks_bl -n_Ys = len(Ys_bl) -n_mus = len(mus_bl) - -# ============================================================================= -# model with uncertain parameters regarding gas extraction -# ============================================================================= -get_uniform_w_frac = lambda b, frac: shape.Uniform(lower=b*(1-frac), upper=b*(1+frac)) - -fset_recir = lambda r: [r, 1] - -def add_degas_params(model, bioreactors, membranes, - recirculation=(1, 19), # ratio to influent Q - b_ermv=0.85, bounds_ermv=(0.5, 1)): - param = model.parameter - # H2E, CH4E = bioreactors - # DM1, DM2 = membranes - lb, ub = recirculation - - for u in bioreactors: - b = 1 - D = shape.Uniform(lb, ub) # equivalent to split varied in [0.5, 0.95] - param(setter=AttrFuncSetter(u, 'split', fset_recir), - name=f'{u.ID}_recirculation_rate', element=u, kind='coupled', - units='', baseline=b, distribution=D) - - # b = b_ermv - # D = shape.Uniform(*bounds_ermv) - # @param(name='H2_removal_efficiency', element=membranes[0], kind='coupled', units='', - # baseline=b, distribution=D) - # def H2_e_rmv(e): - # for dm in membranes: - # dm.H2_degas_efficiency = e - - # D = shape.Uniform(*bounds_ermv) - # @param(name='CH4_removal_efficiency', element=membranes[0], kind='coupled', units='', - # baseline=b, distribution=D) - # def CH4_e_rmv(e): - # for dm in membranes: - # dm.CH4_degas_efficiency = e - - # D = shape.Uniform(*bounds_ermv) - # @param(name='CO2_removal_efficiency', element=DM1, kind='coupled', units='', - # baseline=b, distribution=D) - # def CO2_e_rmv(e): - # DM1.CO2_degas_efficiency = e - # DM2.CO2_degas_efficiency = e - -def add_metrics(model, biogas, wastewater, units): - metric = model.metric - inf, eff = wastewater - R1, R2 = units - S_h2_i_mass = eff.components.S_h2.i_mass - S_ch4_i_mass = eff.components.S_ch4.i_mass - cmps_i_COD = eff.components.i_COD - - # @metric(name='R1 H2', units='mg/L', element='Stage_1') - # def get_H2(): - # return R1.outs[1].iconc['S_h2']*S_h2_i_mass - - # @metric(name='R1 minimum H2 inhibition factor', units='', element='Stage_1') - # def get_Ih2(): - # return min(R1._tempstate['Ih2']) - - # @metric(name='R1 pH', units='', element='Stage_1') - # def get_pH(): - # return R1._tempstate['pH'] - - # @metric(name='R1 pH inhibition factor', units='', element='Stage_1') - # def get_Iph(): - # return R1._tempstate['Iph'][0] - - # @metric(name='R2 CH4', units='mg/L', element='Stage_2') - # def get_CH4(): - # return eff.iconc['S_ch4']*S_ch4_i_mass - - # @metric(name='H2 production', units='kg/d', element='Biogas') - # def get_QH2(): - # return sum([bg.imass['S_h2'] for bg in biogas])*S_h2_i_mass*24 - - # @metric(name='CH4 production', units='kg/d', element='Biogas') - # def get_QCH4(): - # return sum([bg.imass['S_ch4'] for bg in biogas])*S_ch4_i_mass*24 - - @metric(name='Stage 1 COD removal', units='%', element='Stage 1') - def get_rCOD1(): - return (1 - sum(R2.ins[0].mass*cmps_i_COD)/sum(inf.mass*cmps_i_COD))*100 - - @metric(name='Total COD removal', units='%', element='System') - def get_rCOD(): - return (1 - sum(eff.mass*cmps_i_COD)/sum(inf.mass*cmps_i_COD))*100 - - @metric(name='H2 production', units='kgCOD/d', element='Biogas') - def get_QH2(): - return sum([bg.imass['S_h2'] for bg in biogas])*24 - - @metric(name='CH4 production', units='kgCOD/d', element='Biogas') - def get_QCH4(): - return sum([bg.imass['S_ch4'] for bg in biogas])*24 - - # @metric(name='R1 VFAs', units='g/L', element='Stage_1') - # def get_stage1_VFAs(): - # return R1.outs[1].composite('COD', subgroup=vfa_IDs)/1000 - - -def create_modelA(sys=None): - sysA = sys or s.create_systems(which='A')[0] - model = qs.Model(system=sysA, exception_hook='raise') - ws_reg = sysA.flowsheet.stream - inf, eff, bg1, bg2 = ws_reg.BreweryWW_A, ws_reg.Effluent_A, ws_reg.biogas_1A, ws_reg.biogas_2A - u_reg = sysA.flowsheet.unit - H2E, DM1, CH4E, DM2 = u_reg.H2E, u_reg.DM1, u_reg.CH4E, u_reg.DM2 - add_degas_params(model, (H2E, CH4E), (DM1, DM2)) - # add_degas_params(model, (H2E,), (DM1, DM2)) - add_metrics(model, (bg1, bg2), (inf, eff), (H2E, CH4E)) - return model - -def create_modelC(sys=None): - sysC = sys or s.create_systems(which='C')[0] - model = qs.Model(system=sysC, exception_hook='raise') - ws_reg = sysC.flowsheet.stream - inf, eff, bgm1, bgm2, bgh1, bgh2 = ( - ws_reg.BreweryWW_C, - ws_reg.Effluent_C, - ws_reg.biogas_mem_1, - ws_reg.biogas_mem_2, - ws_reg.biogas_hsp_1, - ws_reg.biogas_hsp_2 - ) - u_reg = sysC.flowsheet.unit - R1, DM1, R2, DM2 = u_reg.R1, u_reg.DM1_c, u_reg.R2, u_reg.DM2_c - add_degas_params(model, (R1, R2), (DM1, DM2)) - # add_degas_params(model, (R1,), (DM1, DM2)) - add_metrics(model, (bgm1, bgm2, bgh1, bgh2), (inf, eff), (R1, R2)) - return model - -def create_modelD(sys=None): - sysD = sys or s.create_systems(which='D')[0] - model = qs.Model(system=sysD, exception_hook='raise') - u_reg = sysD.flowsheet.unit - R1, DM1, R2, DM2 = u_reg.R1d, u_reg.DM1d, u_reg.R2d, u_reg.DM2d - ws_reg = sysD.flowsheet.stream - inf, eff, bgm1, bgm2, bgh1, bgh2 = ( - ws_reg.BreweryWW_D, - ws_reg.Effluent_D, - ws_reg.biogas_mem_1d, - ws_reg.biogas_mem_2d, - ws_reg.biogas_hsp_1d, - ws_reg.biogas_hsp_2d - ) - add_degas_params(model, (R1,), (DM1, DM2), recirculation=(1,500)) - add_metrics(model, (bgm1, bgm2, bgh1, bgh2), (inf, eff), (R1, R2)) - return model - -#%% -# ============================================================================= -# model with uncertain ADM1 parameters -# ============================================================================= - -def add_adm_params(model, adm1, units): - param = model.parameter - AnR1, AnR2 = units - for k,v in Ys_bl.items(): - b = v - D = get_uniform_w_frac(b, 0.5) - @param(name=k, element=AnR1, kind='coupled', units='', - baseline=b, distribution=D) - def Y_setter(Y): - pass - - for i in range(len(mus_bl)): - b = mus_bl[i] - D = get_uniform_w_frac(b, 0.5) - @param(name=f'mu_{adm1.IDs[i]}', element=AnR1, kind='coupled', units='d^(-1)', - baseline=b, distribution=D) - def mu_setter(mu): - pass - - for j in range(len(Ks_bl)): - b = Ks_bl[j] - D = get_uniform_w_frac(b, 0.5) - @param(name=f'K_{adm1.IDs[j+4]}', element=AnR1, kind='coupled', units='kg/m3', - baseline=b, distribution=D) - def K_setter(K): - pass - -def create_modelB(sys=None): - sysB = sys or s.create_systems(which='B') - model = qs.Model(system=sysB, exception_hook='raise') - ws_reg = sysB.flowsheet.stream - inf, eff, bg1, bg2 = ws_reg.BreweryWW_B, ws_reg.Effluent_B, ws_reg.biogas_1B, ws_reg.biogas_2B - u_reg = sysB.flowsheet.unit - AnR1, AnR2 = u_reg.AnR1, u_reg.AnR2 - adm1 = AnR1.model - add_adm_params(model, adm1, (AnR1, AnR2)) - add_metrics(model, (bg1, bg2), (inf, eff), (AnR1, AnR2)) - return model - -#%% -# ============================================================================= -# model with random initial conditions -# ============================================================================= -def create_ss_model(system, R1_ID, R2_ID, - R1_baseline_init_conds=None, - R2_baseline_init_conds=None, - frac_var=0.5): - model = qs.Model(system, exception_hook='raise') - _ic1 = R1_baseline_init_conds or s.R1_ss_conds - _ic2 = R2_baseline_init_conds or s.R2_ss_conds - u_reg = system.flowsheet.unit - R1 = getattr(u_reg, R1_ID) - R2 = getattr(u_reg, R2_ID) - param = model.parameter - for ic, u in zip((_ic1,_ic2), (R1,R2)): - for k, v in ic.items(): - b = v - D = get_uniform_w_frac(b, frac_var) - @param(name=f'{k}_0', element=u, kind='coupled', units='mg/L', - baseline=b, distribution=D) - def ic_setter(conc): pass - bgs = [s for s in system.products if s.phase =='g'] - inf, = system.feeds - eff, = set(system.products) - set(bgs) - add_metrics(model, bgs, (inf, eff), (R1, R2)) - return model - -@time_printer -def run_ss_model(model, N, T, t_step, method='BDF', sys_ID=None, - R1_baseline_init_conds=None, - R2_baseline_init_conds=None, - R1_ID='R1', R2_ID='R2', - metrics_path='', timeseries_path='', - rule='L', seed=None, pickle=False): - _ic1 = R1_baseline_init_conds or s.R1_ss_conds - _ic2 = R2_baseline_init_conds or s.R2_ss_conds - k1 = _ic1.keys() - k2 = _ic2.keys() - n_ic1 = len(k1) - u_reg = model._system.flowsheet.unit - R1 = getattr(u_reg, R1_ID) - R2 = getattr(u_reg, R2_ID) - if seed: np.random.seed(seed) - samples = model.sample(N=N, rule=rule) - model.load_samples(samples) - t_span = (0, T) - t_eval = np.arange(0, T+t_step, t_step) - mpath = metrics_path or ospath.join(results_path, f'ss{sys_ID}_table_{seed}.xlsx') - if timeseries_path: tpath = timeseries_path - else: - folder = ospath.join(results_path, f'ss{sys_ID}_time_series_data_{seed}') - os.mkdir(folder) - tpath = ospath.join(folder, 'state.npy') - index = model._index - values = [None] * len(index) - for i, smp in enumerate(samples): - ic1 = dict(zip(k1, smp[:n_ic1])) - ic2 = dict(zip(k2, smp[n_ic1:])) - R1.set_init_conc(**ic1) - R2.set_init_conc(**ic2) - model._system.simulate( - state_reset_hook='reset_cache', - t_span=t_span, - t_eval=t_eval, - method=method, - export_state_to=tpath, - sample_id=i, - ) - values[i] = [i() for i in model.metrics] - model.table[var_indices(model._metrics)] = values - model.table.to_excel(mpath) - R1.set_init_conc(**_ic1) - R2.set_init_conc(**_ic2) - -#%% -@time_printer -def run_model(model, N, T, t_step, method='BDF', - metrics_path='', timeseries_path='', - rule='L', seed=None, pickle=False): - if seed: np.random.seed(seed) - sys_ID = model._system.flowsheet.ID[-1] - samples = model.sample(N=N, rule=rule) - model.load_samples(samples) - t_span = (0, T) - t_eval = np.arange(0, T+t_step, t_step) - mpath = metrics_path or ospath.join(results_path, f'sys{sys_ID}_table_{seed}.xlsx') - if timeseries_path: tpath = timeseries_path - else: - folder = ospath.join(results_path, f'sys{sys_ID}_time_series_data_{seed}') - os.mkdir(folder) - tpath = ospath.join(folder, 'state.npy') - model.evaluate( - state_reset_hook='reset_cache', - t_span=t_span, - t_eval=t_eval, - method=method, - export_state_to=tpath - ) - model.table.to_excel(mpath) - -@time_printer -def run_modelB(model, N, T, t_step, method='BDF', - metrics_path='', timeseries_path='', - rule='L', seed=None, pickle=False): - if seed: np.random.seed(seed) - samples = model.sample(N=N, rule=rule) - model.load_samples(samples) - t_span = (0, T) - t_eval = np.arange(0, T+t_step, t_step) - mpath = metrics_path or ospath.join(results_path, f'sysB_table_{seed}.xlsx') - if timeseries_path: tpath = timeseries_path - else: - folder = ospath.join(results_path, f'sysB_time_series_data_{seed}') - os.mkdir(folder) - tpath = ospath.join(folder, 'state.npy') - index = model._index - values = [None] * len(index) - adm1 = model._system.flowsheet.unit.AnR1.model - for i in index: - smp = samples[i] - Ys = dict(zip(Ys_bl.keys(), smp[:n_Ys])) - adm1.set_parameters(**Ys) - adm1.rate_function._params['rate_constants'][:] = smp[n_Ys: (n_Ys+n_mus)] - adm1.rate_function._params['half_sat_coeffs'][:] = smp[(n_Ys+n_mus):] - model._system.simulate( - state_reset_hook='reset_cache', - t_span=t_span, - t_eval=t_eval, - method=method, - export_state_to=tpath, - sample_id=i, - ) - values[i] = [i() for i in model.metrics] - model.table[var_indices(model._metrics)] = values - model.table.to_excel(mpath) \ No newline at end of file diff --git a/exposan/metab_mock/graveyard/systems.py b/exposan/metab_mock/graveyard/systems.py deleted file mode 100644 index 07f6b13e..00000000 --- a/exposan/metab_mock/graveyard/systems.py +++ /dev/null @@ -1,450 +0,0 @@ -# -*- coding: utf-8 -*- -''' -EXPOsan: Exposition of sanitation and resource recovery systems - -This module is developed by: - - Joy Zhang - -This module is under the University of Illinois/NCSA Open Source License. -Please refer to https://github.com/QSD-Group/EXPOsan/blob/main/LICENSE.txt -for license details. -''' - -import numpy as np -import qsdsan as qs -from qsdsan import sanunits as su, processes as pc, WasteStream, System -from qsdsan.utils import time_printer, ospath -from chemicals.elements import molecular_weight as get_mw -from exposan.metab_mock import DegassingMembrane as DM, METAB_AnCSTR as AB - -folder = ospath.dirname(__file__) - -__all__ = ( - 'create_systems', - 'default_inf_concs', - 'default_R1_init_conds', - 'default_R2_init_conds', - 'R1_ss_conds', - 'R2_ss_conds', - 'yields_bl', 'mus_bl', 'Ks_bl', - 'biomass_IDs', - 'vfa_IDs' - ) - -#%% default values -Q = 5 # influent flowrate [m3/d] -T1 = 273.15+35 # temperature [K] -Vl1 = 5 # liquid volume [m^3] -Vg1 = 0.556 # headspace volume [m^3] -split_1 = 0.75 # split ratio to side-stream -tau_1 = 0.021 # degassing membrane retention time [d] - -T2 = 273.15+25 -Vl2 = 75 -Vg2 = 5 -split_2 = 0.75 -tau_2 = 0.021 - -fermenters = ('X_su', 'X_aa', 'X_fa', 'X_c4', 'X_pro') -methanogens = ('X_ac', 'X_h2') -biomass_IDs = (*fermenters, *methanogens) -vfa_IDs = ('S_va', 'S_bu', 'S_pro', 'S_ac') - -C_mw = get_mw({'C':1}) -N_mw = get_mw({'N':1}) - -default_inf_concs = { - 'S_su':3.0, - 'S_aa':0.6, - 'S_fa':0.4, - 'S_va':0.4, - 'S_bu':0.4, - 'S_pro':0.4, - 'S_ac':0.4, - 'S_h2':5e-9, - 'S_ch4':5e-6, - 'S_IC':0.04*C_mw, - 'S_IN':0.01*N_mw, - 'S_I':0.02, - 'X_c':0.1, - 'X_ch':0.3, - 'X_pr':0.5, - 'X_li':0.25, - 'X_aa':1e-3, - 'X_fa':1e-3, - 'X_c4':1e-3, - 'X_pro':1e-3, - 'X_ac':1e-3, - 'X_h2':1e-3, - 'X_I':0.025, - 'S_cat':0.04, - 'S_an':0.02 - } - -yields_bl = { - 'Y_su': 0.1, - 'Y_aa': 0.08, - 'Y_fa': 0.06, - 'Y_c4': 0.06, - 'Y_pro': 0.04, - 'Y_ac': 0.05, - 'Y_h2': 0.06 - } - -mus_bl = np.array([5.0e-01, 1.0e+01, 1.0e+01, 1.0e+01, 3.0e+01, 5.0e+01, 6.0e+00, - 2.0e+01, 2.0e+01, 1.3e+01, 8.0e+00, 3.5e+01, 2.0e-02, 2.0e-02, - 2.0e-02, 2.0e-02, 2.0e-02, 2.0e-02, 2.0e-02]) - -Ks_bl = np.array([5.0e-01, 3.0e-01, 4.0e-01, 2.0e-01, - 2.0e-01, 1.0e-01, 1.5e-01, 7.0e-06]) - -default_R1_init_conds = { - 'S_su': 0.0124*1e3, - 'S_aa': 0.0055*1e3, - 'S_fa': 0.1074*1e3, - 'S_va': 0.0123*1e3, - 'S_bu': 0.0140*1e3, - 'S_pro': 0.0176*1e3, - 'S_ac': 0.0893*1e3, - 'S_h2': 2.5055e-7*1e3, - 'S_ch4': 0.0555*1e3, - 'S_IC': 0.0951*C_mw*1e3, - 'S_IN': 0.0945*N_mw*1e3, - 'S_I': 0.1309*1e3, - 'X_ch': 0.0205*1e3, - 'X_pr': 0.0842*1e3, - 'X_li': 0.0436*1e3, - 'X_su': 1.87*1e3, - 'X_aa': 5.58*1e3, - 'X_fa': 2.03*1e3, - 'X_c4': 2.15*1e3, - 'X_pro': 1.00*1e3, - } - -default_R2_init_conds = { - 'S_su': 0.0124*1e3, - 'S_aa': 0.0055*1e3, - 'S_fa': 0.1074*1e3, - 'S_va': 0.0123*1e3, - 'S_bu': 0.0140*1e3, - 'S_pro': 0.0176*1e3, - 'S_ac': 0.0893*1e3, - 'S_h2': 2.5055e-7*1e3, - 'S_ch4': 0.0555*1e3, - 'S_IC': 0.0951*C_mw*1e3, - 'S_IN': 0.0945*N_mw*1e3, - 'S_I': 0.1309*1e3, - 'X_ac': 8.80*1e3, - 'X_h2': 3.70*1e3, - } - -R1_ss_conds = { - 'S_su': 0.0145871088552909*1e3, - 'S_aa': 0.00643308564144693*1e3, - 'S_fa': 0.634823005990967*1e3, - 'S_va': 0.624510322247682*1e3, - 'S_bu': 1.03793927591996*1e3, - 'S_pro': 1.24676871525373*1e3, - 'S_ac': 2.00250371674824*1e3, - 'S_h2': 0.00850364943532684*1e3, - 'S_ch4': 0.0000422133982597226*1e3, - 'S_IC': 0.0951*C_mw*1e3, - 'S_IN': 0.0945*N_mw*1e3, - 'S_I': 0.027310256066728*1e3, - 'X_c': 0.146203507058736*1e3, - 'X_ch': 0.0286018513139117*1e3, - 'X_pr': 0.0467836694957302*1e3, - 'X_li': 0.0247209587890493*1e3, - 'X_su': 4.69052782535406*1e3, - 'X_aa': 1.22829926704024*1e3, - 'X_fa': 0.0147446263753011*1e3, - 'X_c4': 0.0149933579422897*1e3, - 'X_pro': 0.0145343147735253*1e3, - 'X_ac': 0.00098041337766024*1e3, - 'X_h2': 0.00110808891184369*1e3, - 'X_I': 0.0396205121367899*1e3 - } - -R2_ss_conds = { - 'S_su': 0.00106990968535691*1e3, - 'S_aa': 0.00125571416517827*1e3, - 'S_fa': 0.121097573221394*1e3, - 'S_va': 0.0132519103137696*1e3, - 'S_bu': 0.0172912281196732*1e3, - 'S_pro': 0.020032163173878*1e3, - 'S_ac': 0.00574002755366853*1e3, - 'S_h2': 3.76969944940856e-08*1e3, - 'S_ch4': 0.0499411746585487*1e3, - 'S_IC': 0.0951*C_mw*1e3, - 'S_IN': 0.0945*N_mw*1e3, - 'S_I': 0.105601391746794*1e3, - 'X_c': 0.0897520281015078*1e3, - 'X_ch': 0.00108163641708242*1e3, - 'X_pr': 0.00120204580901502*1e3, - 'X_li': 0.00150204523369107*1e3, - 'X_su': 0.195961987850137*1e3, - 'X_aa': 0.059723477130333*1e3, - 'X_fa': 0.0351858744892462*1e3, - 'X_c4': 0.0812315951844566*1e3, - 'X_pro': 0.0503466475437059*1e3, - 'X_ac': 1.1653549028287*1e3, - 'X_h2': 0.4352809013846*1e3, - 'X_I': 0.196117291164614*1e3 - } - -# # R1_split = (0.9, 0.1), <0.2% diff if uses these -# R1_ss_conds = { -# 'S_su': 0.01431332653033124*1e3, -# 'S_aa': 0.006319229324446604*1e3, -# 'S_fa': 0.6348272949518201*1e3, -# 'S_va': 0.6245549045932813*1e3, -# 'S_bu': 1.038036275314748*1e3, -# 'S_pro': 1.2468433077261327*1e3, -# 'S_ac': 2.0025986879937308*1e3, -# 'S_h2': 0.011117814005879386*1e3, -# 'S_ch4': 5.952203539352817e-05*1e3, -# 'S_IC': 0.08330274953356481*C_mw*1e3, -# 'S_IN': 0.21612074810295764*N_mw*1e3, -# 'S_I': 0.02730942113623604*1e3, -# 'X_c': 0.14618842266689933*1e3, -# 'X_ch': 0.028601712932738325*1e3, -# 'X_pr': 0.04678353111459632*1e3, -# 'X_li': 0.024720751217894823*1e3, -# 'X_su': 4.6915202551212944*1e3, -# 'X_aa': 1.2274469593168187*1e3, -# 'X_fa': 0.014302724094919613*1e3, -# 'X_c4': 0.014436386133097296*1e3, -# 'X_pro': 0.014311685377540902*1e3, -# 'X_ac': 0.0009804172325111186*1e3, -# 'X_h2': 0.0011332727803773114*1e3, -# 'X_I': 0.039618842271433245*1e3 -# } - -# R2_ss_conds = { -# 'S_su': 0.0008282629572271559*1e3, -# 'S_aa':0.001032694581229776*1e3, -# 'S_fa': 0.12736791568063396*1e3, -# 'S_va': 0.013392575245557583*1e3, -# 'S_bu': 0.01747268169433745*1e3, -# 'S_pro': 0.020406491318645616*1e3, -# 'S_ac': 0.00999472311030857*1e3, -# 'S_h2': 7.706131429350044e-08*1e3, -# 'S_ch4': 0.07044765251475965*1e3, -# 'S_IC': 0.49869691792599646*C_mw*1e3, -# 'S_IN': 0.20439321199252136*N_mw*1e3, -# 'S_I': 0.07350288752133147*1e3, -# 'X_c': 0.06137141129938335*1e3, -# 'X_ch': 0.0007989050950664583*1e3, -# 'X_pr': 0.0009192059152983118*1e3, -# 'X_li': 0.0010780568238771538*1e3, -# 'X_su': 0.19130327796202837*1e3, -# 'X_aa': 0.056026151032630665*1e3, -# 'X_fa': 0.03107211586931466*1e3, -# 'X_c4': 0.07959903792107152*1e3, -# 'X_pro': 0.04955864019101192*1e3, -# 'X_ac': 0.6361761447275743*1e3, -# 'X_h2':0.20772282044789275*1e3, -# 'X_I': 0.13200577400466326*1e3 -# } - - -# %% -# ============================================================================= -# Preliminary analyses with mock METAB configuration -# ============================================================================= - -def create_systems(flowsheet_A=None, flowsheet_B=None, flowsheet_C=None, - flowsheet_D=None, flowsheet_E=None, - inf_concs={}, R1_init_conds={}, R2_init_conds={}, which=None): - which = which or ('A', 'B', 'C', 'D', 'E') - if isinstance(which, str): which = (which,) - ############# load components and set thermo ############# - pc.create_adm1_cmps() - inf_concs = inf_concs or default_inf_concs - brewery_ww = WasteStream('BreweryWW_A', T=T1) - brewery_ww.set_flow_by_concentration(Q, concentrations=inf_concs, units=('m3/d', 'kg/m3')) - ############# load process model ########################### - adm1 = pc.ADM1() - R1_init_conds = R1_init_conds or default_R1_init_conds - R2_init_conds = R2_init_conds or default_R2_init_conds - systems = [] - - if 'A' in which: - flowsheet_A = flowsheet_A or qs.Flowsheet('METAB_sysA') - qs.main_flowsheet.set_flowsheet(flowsheet_A) - ############# create WasteStream objects ################# - eff_A = WasteStream('Effluent_A', T=T2) - bg1_A = WasteStream('biogas_1A', phase='g') - bg2_A = WasteStream('biogas_2A', phase='g') - flowsheet_A.stream.register(brewery_ww.ID, brewery_ww) - - ############# sysA unit operation ######################## - H2E = AB('H2E', ins=[brewery_ww, 'return_1'], outs=('sidestream_1', ''), - split=(split_1, 1-split_1), V_liq=Vl1, V_gas=Vg1, T=T1, model=adm1, - retain_cmps=fermenters) - DM1 = DM('DM1', ins=H2E-0, outs=(bg1_A, 1-H2E), tau=tau_1) - CH4E = AB('CH4E', ins=[H2E-1, 'return_2'], outs=('sidestream_2', eff_A), - split=(split_2, 1-split_2), V_liq=Vl2, V_gas=Vg2, T=T2, model=adm1, - retain_cmps=methanogens) - DM2 = DM('DM2', ins=CH4E-0, outs=(bg2_A, 1-CH4E), tau=tau_2) - H2E.set_init_conc(**R1_ss_conds) - CH4E.set_init_conc(**R2_ss_conds) - # H2E.set_init_conc(**R1_init_conds) - # CH4E.set_init_conc(**R2_init_conds) - sysA = System('mock_METAB', - path=(H2E, DM1, CH4E, DM2), - recycle=(DM1-1, DM2-1)) - sysA.set_dynamic_tracker(H2E, CH4E, bg1_A, bg2_A) - systems.append(sysA) - - if 'B' in which: - flowsheet_B = flowsheet_B or qs.Flowsheet('METAB_sysB') - qs.main_flowsheet.set_flowsheet(flowsheet_B) - - ############# sysB streams ######################## - inf_b = brewery_ww.copy('BreweryWW_B') - eff_B = WasteStream('Effluent_B', T=T2) - bg1_B = WasteStream('biogas_1B', phase='g') - bg2_B = WasteStream('biogas_2B', phase='g') - - ############# sysB unit operation ################# - AnR1 = su.AnaerobicCSTR('AnR1', ins=inf_b, outs=(bg1_B, ''), - V_liq=Vl1, V_gas=Vg1, T=T1, model=adm1, - retain_cmps=fermenters) - AnR2 = su.AnaerobicCSTR('AnR2', ins=AnR1-1, outs=(bg2_B, eff_B), - V_liq=Vl2, V_gas=Vg2, T=T2, model=adm1, - retain_cmps=methanogens) - # AnR1.set_init_conc(**R1_init_conds) - # AnR2.set_init_conc(**R2_init_conds) - AnR1.set_init_conc(**R1_ss_conds) - AnR2.set_init_conc(**R2_ss_conds) - sysB = System('baseline', path=(AnR1, AnR2)) - sysB.set_dynamic_tracker(AnR1, AnR2, bg1_B, bg2_B) - systems.append(sysB) - - if 'C' in which: - flowsheet_C = flowsheet_C or qs.Flowsheet('METAB_sysC') - qs.main_flowsheet.set_flowsheet(flowsheet_C) - - ############# sysC streams ######################## - inf_c = brewery_ww.copy('BreweryWW_C') - eff_c = WasteStream('Effluent_C', T=T2) - bgm1 = WasteStream('biogas_mem_1', phase='g') - bgm2 = WasteStream('biogas_mem_2', phase='g') - bgh1 = WasteStream('biogas_hsp_1', phase='g') - bgh2 = WasteStream('biogas_hsp_2', phase='g') - - ############# sysC unit operation ################# - R1 = su.AnaerobicCSTR('R1', ins=[inf_c, 'return_1'], - outs=(bgh1, 'sidestream_1', ''), - split=(split_1, 1-split_1), - V_liq=Vl1, V_gas=Vg1, T=T1, model=adm1, - retain_cmps=fermenters) - DM1c = DM('DM1_c', ins=R1-1, outs=(bgm1, 1-R1), tau=tau_1) - # DM1c = DM('DM1_c', ins=R1-1, outs=(bgm1, 1-R1), tau=0.1) - - R2 = su.AnaerobicCSTR('R2', ins=[R1-2, 'return_2'], - outs=(bgh2, 'sidestream_2', eff_c), - split=(split_2, 1-split_2), - V_liq=Vl2, V_gas=Vg2, T=T2, model=adm1, - retain_cmps=methanogens) - DM2c = DM('DM2_c', ins=R2-1, outs=(bgm2, 1-R2), tau=tau_2) - # DM2c = DM('DM2_c', ins=R2-1, outs=(bgm2, 1-R2), tau=0.1) - R1.set_init_conc(**R1_ss_conds) - R2.set_init_conc(**R2_ss_conds) - sysC = System('combined_METAB', path=(R1, DM1c, R2, DM2c), - recycle=(DM1c-1, DM2c-1)) - sysC.set_dynamic_tracker(R1, R2, bgm1, bgm2, bgh1, bgh2) - systems.append(sysC) - - if 'D' in which: - flowsheet_D = flowsheet_D or qs.Flowsheet('METAB_sysD') - qs.main_flowsheet.set_flowsheet(flowsheet_D) - - ############# sysC streams ######################## - inf_d = brewery_ww.copy('BreweryWW_D') - eff_d = WasteStream('Effluent_D', T=T2) - bgm1_d = WasteStream('biogas_mem_1d', phase='g') - bgm2_d = WasteStream('biogas_mem_2d', phase='g') - bgh1_d = WasteStream('biogas_hsp_1d', phase='g') - bgh2_d = WasteStream('biogas_hsp_2d', phase='g') - - ############# sysC unit operation ################# - R1d = su.AnaerobicCSTR('R1d', ins=[inf_d, 'return_1d'], - outs=(bgh1_d, 'sidestream_1d', ''), - split=(split_1, 1-split_1), - V_liq=Vl1, V_gas=Vg1, T=T1, model=adm1, - retain_cmps=fermenters) - DM1d = DM('DM1d', ins=R1d-1, outs=(bgm1_d, 1-R1d), tau=tau_1) - - R2d = su.AnaerobicCSTR('R2d', ins=R1d-2, - outs=(bgh2_d, ''), - V_liq=Vl2, V_gas=Vg2, T=T2, model=adm1, - retain_cmps=methanogens) - DM2d = DM('DM2d', ins=R2d-1, outs=(bgm2_d, eff_d), tau=tau_2) - R1d.set_init_conc(**R1_ss_conds) - R2d.set_init_conc(**R2_ss_conds) - sysD = System('hybrid_METAB', path=(R1d, DM1d, R2d, DM2d), - recycle=(DM1d-1, )) - sysD.set_dynamic_tracker(R1d, R2d, bgm1_d, bgm2_d, bgh1_d, bgh2_d) - systems.append(sysD) - - if 'E' in which: - flowsheet_E = flowsheet_E or qs.Flowsheet('METAB_sysE') - qs.main_flowsheet.set_flowsheet(flowsheet_E) - - ############# sysC streams ######################## - inf_e = brewery_ww.copy('BreweryWW_E') - eff_e = WasteStream('Effluent_E', T=T2) - # bgm1_d = WasteStream('biogas_mem_1d', phase='g') - bgm2_e = WasteStream('biogas_mem_2e', phase='g') - bgh1_e = WasteStream('biogas_hsp_1e', phase='g') - bgh2_e = WasteStream('biogas_hsp_2e', phase='g') - - ############# sysC unit operation ################# - R1e = su.AnaerobicCSTR('R1e', ins=inf_e, outs=(bgh1_e, ''), - fixed_headspace_P=True, - headspace_P=0.1013, - V_liq=Vl1, V_gas=Vg1, T=T1, model=adm1, - retain_cmps=fermenters) - R2e = su.AnaerobicCSTR('R2e', ins=R1e-1, outs=(bgh2_e, ''), - V_liq=Vl2, V_gas=Vg2, T=T2, model=adm1, - retain_cmps=methanogens) - DM2e = DM('DM2e', ins=R2e-1, outs=(bgm2_e, eff_e), tau=tau_2) - R1e.set_init_conc(**R1_ss_conds) - R2e.set_init_conc(**R2_ss_conds) - sysE = System('Best_METAB', path=(R1e, R2e, DM2e)) - sysE.set_dynamic_tracker(R1e, R2e, bgm2_e, bgh1_e, bgh2_e) - systems.append(sysE) - - return systems - -#%% -@time_printer -def run(t, t_step, method=None, **kwargs): - global sysA, sysB, sysC, sysD, sysE - sysA, sysB, sysC, sysD, sysE = create_systems() - for sys in (sysA, sysB, sysC, sysD): - print(f'Simulating {sys.ID}...') - sys.simulate(state_reset_hook='reset_cache', - t_span=(0,t), - t_eval=np.arange(0, t+t_step, t_step), - method=method, - # export_state_to=ospath.join(folder, f'results/{method}_{t}d_{sys.ID[-4:]}.xlsx'), - **kwargs) - - -if __name__ == '__main__': - t = 200 - t_step = 5 - # method = 'RK45' - # method = 'RK23' - # method = 'DOP853' - # method = 'Radau' - method = 'BDF' - # method = 'LSODA' - msg = f'Method {method}' - print(f'\n{msg}\n{"-"*len(msg)}') # long live OCD! - print(f'Time span 0-{t}d \n') - run(t, t_step, method=method) \ No newline at end of file diff --git a/exposan/metab_mock/hrt_test.py b/exposan/metab_mock/hrt_test.py deleted file mode 100644 index 43df3a33..00000000 --- a/exposan/metab_mock/hrt_test.py +++ /dev/null @@ -1,96 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Wed Jan 25 12:48:54 2023 - -@author: joy_c -""" - -from qsdsan import processes as pc, sanunits as su, WasteStream, System -from qsdsan.utils import ExogenousDynamicVariable as EDV -from exposan.metab_mock import ( - rhos_adm1_ph_ctrl, - default_inf_concs, - default_R1_init_conds, - R1_ss_conds, - biomass_IDs, - ) -import numpy as np - -pc.create_adm1_cmps() -adm1 = pc.ADM1() -# dyn_params = adm1.rate_function.params.copy() -# adm1.set_rate_function(rhos_adm1_ph_ctrl) -# adm1.rate_function._params = dyn_params - -Q = 5 -T1 = 298.15 -T2 = 298.15 -ph1 = 5.8 -ph2 = 7.2 - -Temp = EDV('T1', function=lambda t: T1) -Temp2 = EDV('T2', function=lambda t: T2) -pH1 = EDV('pH1', function=lambda t: ph1) -pH2 = EDV('pH2', function=lambda t: ph2) - -ics = { - 'S_su': 0.32679646003805314, - 'S_aa': 0.3801610819236484, - 'S_fa': 12.437222319633748, - 'S_va': 0.3719673543175543, - 'S_bu': 0.47283246583627593, - 'S_pro': 0.3946420365926535, - 'S_ac': 10.182894473261367, - 'S_h2': 1.1655700001506622e-05, - 'S_ch4': 67.17348627201263, - 'S_IC': 846.4879614661522, - 'S_IN': 222.13725282096297, - 'S_I': 110.71467278942289, - 'X_c': 107.43132381172228, - 'X_ch': 1.2600235711799973, - 'X_pr': 1.3804329631122664, - 'X_li': 1.7696259648387357, - 'X_su': 732.9760678333023, - 'X_aa': 224.81751931525334, - 'X_fa': 126.7301174776879, - 'X_c4': 227.8726398428066, - 'X_pro': 140.2738127019708, - 'X_ac': 669.4626559278454, - 'X_h2': 245.67774602566578, - 'X_I': 206.42934561053158, - 'S_cat': 40.0, - 'S_an': 20.0, - } - -inf = WasteStream('inf') -eff = WasteStream('eff') -bg = WasteStream('bg', phase='g') -inf.set_flow_by_concentration(Q, concentrations=default_inf_concs, - units=('m3/d', 'kg/m3')) - -R1 = su.AnaerobicCSTR('R1', ins=inf, outs=(bg, eff), - # V_liq=50, V_gas=5, - V_liq=20, V_gas=2, - model=adm1, T=T1, - retain_cmps=biomass_IDs, - # exogenous_vars=(Temp1, pH1), - ) -# R1.set_init_conc(**R1_ss_conds) -R1.set_init_conc(**ics) - -sys = System('sys', path=(R1,)) -sys.set_dynamic_tracker(R1, eff, bg) - -#%% -HRTs = 2**np.linspace(-1, 4) -rCOD = {} -for i, tau in enumerate(HRTs): - R1.V_liq = Q*tau - R1.V_gas = R1.V_liq/(i+8) - sys.simulate(state_reset_hook='reset_cache', - t_span=(0, 200), - method='BDF' - ) - assert np.all(R1._state >= 0) - R1.scope.plot_time_series(('S_h2', 'S_ch4')) - rCOD[tau] = 1-eff.COD/inf.COD \ No newline at end of file diff --git a/exposan/metab_mock/init_cond_test.py b/exposan/metab_mock/init_cond_test.py deleted file mode 100644 index bdec7782..00000000 --- a/exposan/metab_mock/init_cond_test.py +++ /dev/null @@ -1,33 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Tue Feb 7 16:55:35 2023 - -@author: joy_c -""" - -from exposan.metab_mock import create_systems, biomass_IDs -import numpy as np - -sys, = create_systems(which='C') -u = sys.units[0] -bg, eff = u.outs -inf, = u.ins -idx = inf.components.indices(biomass_IDs) -_ic = u._concs[idx].copy() -init_bm = np.linspace(1, 10, 10) - -u.V_liq = 2.5 -u.V_gas = 0.25 - -rcod = [] -bm = [] -bm_dist = [] -fug_ch4 = [] - -for i in init_bm: - u._concs[idx] = _ic*i - sys.simulate(state_reset_hook='reset_cache', t_span=(0, 400), method='BDF') - rcod.append(1-eff.COD/inf.COD) - bm.append(sum(u._state[idx])) - bm_dist.append(u._state[idx]/sum(u._state[idx])*100) - fug_ch4.append(eff.iconc['S_ch4']) diff --git a/exposan/metab_mock/process.py b/exposan/metab_mock/process.py deleted file mode 100644 index c6a023ff..00000000 --- a/exposan/metab_mock/process.py +++ /dev/null @@ -1,99 +0,0 @@ -# -*- coding: utf-8 -*- -''' -EXPOsan: Exposition of sanitation and resource recovery systems - -This module is developed by: - - Joy Zhang - -This module is under the University of Illinois/NCSA Open Source License. -Please refer to https://github.com/QSD-Group/EXPOsan/blob/main/LICENSE.txt -for license details. -''' - -import numpy as np -from math import log10 -from scipy.optimize import brenth -from qsdsan.processes._adm1 import ( - R, - mass2mol_conversion, - T_correction_factor, - acid_base_rxn, - substr_inhibit, - Hill_inhibit, - non_compet_inhibit - ) - -__all__ = ('rhos_adm1_ph_ctrl',) - - -#%% -rhos = np.zeros(22) # 22 kinetic processes -Cs = np.empty(19) - -def rhos_adm1_ph_ctrl(state_arr, params): - ks = params['rate_constants'] - Ks = params['half_sat_coeffs'] - cmps = params['components'] - pH_ULs = params['pH_ULs'] - pH_LLs = params['pH_LLs'] - KS_IN = params['KS_IN'] - KI_nh3 = params['KI_nh3'] - KIs_h2 = params['KIs_h2'] - KHb = params['K_H_base'] - Kab = params['Ka_base'] - KH_dH = params['K_H_dH'] - Ka_dH = params['Ka_dH'] - kLa = params['kLa'] - T_base = params['T_base'] - root = params['root'] - - Cs[:8] = state_arr[12:20] - Cs[8:12] = state_arr[19:23] - Cs[12:] = state_arr[16:23] - substrates = state_arr[:8] - S_va, S_bu, S_h2, S_IN = state_arr[[3,4,7,10]] - unit_conversion = mass2mol_conversion(cmps) - cmps_in_M = state_arr[:27] * unit_conversion - weak_acids = cmps_in_M[[24, 25, 10, 9, 6, 5, 4, 3]] - - T_op, pH = state_arr[-2:] #!!! change in state_arr - biogas_S = state_arr[7:10].copy() - biogas_p = R * T_op * state_arr[27:30] - Kas = Kab * T_correction_factor(T_base, T_op, Ka_dH) - KH = KHb * T_correction_factor(T_base, T_op, KH_dH) / unit_conversion[7:10] - - rhos[:-3] = ks * Cs - Monod = substr_inhibit(substrates, Ks) - rhos[4:12] *= Monod - if S_va > 0: rhos[7] *= 1/(1+S_bu/S_va) - if S_bu > 0: rhos[8] *= 1/(1+S_va/S_bu) - - h = 10**(-pH) - delta = acid_base_rxn(h, weak_acids, Kas) - S_cat = weak_acids[0] - delta - - nh3 = Kas[1] * weak_acids[2] / (Kas[1] + h) - co2 = weak_acids[3] - Kas[2] * weak_acids[3] / (Kas[2] + h) - biogas_S[-1] = co2 / unit_conversion[9] - - Iph = Hill_inhibit(h, pH_ULs, pH_LLs) - Iin = substr_inhibit(S_IN, KS_IN) - Ih2 = non_compet_inhibit(S_h2, KIs_h2) - Inh3 = non_compet_inhibit(nh3, KI_nh3) - rhos[4:12] *= Iph * Iin - rhos[6:10] *= Ih2 - rhos[10] *= Inh3 - rhos[-3:] = kLa * (biogas_S - KH * biogas_p) - root.data = { - 'pH':pH, - 'Iph':Iph, - 'Ih2':Ih2, - 'Iin':Iin, - 'Inh3':Inh3, - 'Monod':Monod, - 'rhos':rhos[4:12].copy(), - 'gas_transfer':rhos[-3:].copy(), - 'S_cat':S_cat - } - return rhos \ No newline at end of file diff --git a/exposan/metab_mock/systems.py b/exposan/metab_mock/systems.py deleted file mode 100644 index 929902cc..00000000 --- a/exposan/metab_mock/systems.py +++ /dev/null @@ -1,412 +0,0 @@ -# -*- coding: utf-8 -*- -''' -EXPOsan: Exposition of sanitation and resource recovery systems - -This module is developed by: - - Joy Zhang - -This module is under the University of Illinois/NCSA Open Source License. -Please refer to https://github.com/QSD-Group/EXPOsan/blob/main/LICENSE.txt -for license details. -''' - -import numpy as np, qsdsan as qs -from qsdsan import ( - processes as pc, - WasteStream, System, TEA, LCA, PowerUtility, - ImpactItem as IItm, - StreamImpactItem as SIItm, - ) -from qsdsan.utils import ExogenousDynamicVariable as EDV -from exposan.metab_mock import ( - _impact_item_loaded, - load_lca_data, - rhos_adm1_ph_ctrl, - METAB_AnCSTR as AB - ) -from exposan.metab import ( - DegassingMembrane as DM, - IronSpongeTreatment as IST, - DoubleMembraneGasHolder as GH, - add_strm_iitm, add_TEA_LCA - ) - -__all__ = ( - 'create_systems', - 'default_inf_concs', - 'default_R1_init_conds', - 'default_R2_init_conds', - 'default_R_init_conds', - 'R1_ss_conds', - 'R2_ss_conds', - 'yields_bl', 'mus_bl', 'Ks_bl', - 'fermenters', 'methanogens', 'biomass_IDs', - 'vfa_IDs', - ) - -#%% default values -scale = 1 -Q = 5*scale # influent flowrate [m3/d] -T1 = 273.15+35 # temperature [K] -# T1 = 273.15+25 -Vl1 = 5*scale # liquid volume [m^3] -Vg1 = 0.556*(scale**0.5) # headspace volume [m^3] -ph1 = 5.8 - -T2 = 273.15+25 -Vl2 = 75*scale -Vg2 = 5*(scale**0.5) -ph2 = 7.2 - -T3 = T2 -Vl3 = 5*scale -Vg3 = 0.556*(scale**0.5) - -bl = 1 # yr, bead lifetime -# bl = 10 - -Temp1 = EDV('T1', function=lambda t: T1) -Temp2 = EDV('T2', function=lambda t: T2) -pH1 = EDV('pH1', function=lambda t: ph1) -pH2 = EDV('pH2', function=lambda t: ph2) - -fermenters = ('X_su', 'X_aa', 'X_fa', 'X_c4', 'X_pro') -methanogens = ('X_ac', 'X_h2') -biomass_IDs = (*fermenters, *methanogens) -vfa_IDs = ('S_va', 'S_bu', 'S_pro', 'S_ac') - -C_mw = 12.0107 -N_mw = 14.0067 - -default_inf_concs = { - 'S_su':3.0, - 'S_aa':0.6, - 'S_fa':0.4, - 'S_va':0.4, - 'S_bu':0.4, - 'S_pro':0.4, - 'S_ac':0.4, - 'S_h2':5e-9, - 'S_ch4':5e-6, - 'S_IC':0.04*C_mw, - 'S_IN':0.01*N_mw, - 'S_I':0.02, - 'X_c':0.1, - 'X_ch':0.3, - 'X_pr':0.5, - 'X_li':0.25, - 'X_aa':1e-3, - 'X_fa':1e-3, - 'X_c4':1e-3, - 'X_pro':1e-3, - 'X_ac':1e-3, - 'X_h2':1e-3, - 'X_I':0.025, - 'S_cat':0.04, - 'S_an':0.02 - } - -yields_bl = { - 'Y_su': 0.1, - 'Y_aa': 0.08, - 'Y_fa': 0.06, - 'Y_c4': 0.06, - 'Y_pro': 0.04, - 'Y_ac': 0.05, - 'Y_h2': 0.06 - } - -mus_bl = np.array([5.0e-01, 1.0e+01, 1.0e+01, 1.0e+01, 3.0e+01, 5.0e+01, 6.0e+00, - 2.0e+01, 2.0e+01, 1.3e+01, 8.0e+00, 3.5e+01, 2.0e-02, 2.0e-02, - 2.0e-02, 2.0e-02, 2.0e-02, 2.0e-02, 2.0e-02]) - -Ks_bl = np.array([5.0e-01, 3.0e-01, 4.0e-01, 2.0e-01, - 2.0e-01, 1.0e-01, 1.5e-01, 7.0e-06]) - -default_R1_init_conds = { - 'S_su': 0.0124*1e3, - 'S_aa': 0.0055*1e3, - 'S_fa': 0.1074*1e3, - 'S_va': 0.0123*1e3, - 'S_bu': 0.0140*1e3, - 'S_pro': 0.0176*1e3, - 'S_ac': 0.0893*1e3, - 'S_h2': 2.5055e-7*1e3, - 'S_ch4': 0.0555*1e3, - 'S_IC': 0.0951*C_mw*1e3, - 'S_IN': 0.0945*N_mw*1e3, - 'S_I': 0.1309*1e3, - 'X_ch': 0.0205*1e3, - 'X_pr': 0.0842*1e3, - 'X_li': 0.0436*1e3, - 'X_su': 1.87*1e3, - 'X_aa': 5.58*1e3, - 'X_fa': 2.03*1e3, - 'X_c4': 2.15*1e3, - 'X_pro': 1.00*1e3, - } - -default_R2_init_conds = { - 'S_su': 0.0124*1e3, - 'S_aa': 0.0055*1e3, - 'S_fa': 0.1074*1e3, - 'S_va': 0.0123*1e3, - 'S_bu': 0.0140*1e3, - 'S_pro': 0.0176*1e3, - 'S_ac': 0.0893*1e3, - 'S_h2': 2.5055e-7*1e3, - 'S_ch4': 0.0555*1e3, - 'S_IC': 0.0951*C_mw*1e3, - 'S_IN': 0.0945*N_mw*1e3, - 'S_I': 0.1309*1e3, - 'X_ac': 8.80*1e3, - 'X_h2': 3.70*1e3, - } - -default_R_init_conds = { - 'S_su': 0.32679646003805314, - 'S_aa': 0.3801610819236484, - 'S_fa': 12.437222319633748, - 'S_va': 0.3719673543175543, - 'S_bu': 0.47283246583627593, - 'S_pro': 0.3946420365926535, - 'S_ac': 10.182894473261367, - 'S_h2': 1.1655700001506622e-05, - 'S_ch4': 67.17348627201263, - 'S_IC': 846.4879614661522, - 'S_IN': 222.13725282096297, - 'S_I': 110.71467278942289, - 'X_c': 107.43132381172228, - 'X_ch': 1.2600235711799973, - 'X_pr': 1.3804329631122664, - 'X_li': 1.7696259648387357, - 'X_su': 732.9760678333023*5, - 'X_aa': 224.81751931525334*5, - 'X_fa': 126.7301174776879*5, - 'X_c4': 227.8726398428066*5, - 'X_pro': 140.2738127019708*5, - 'X_ac': 669.4626559278454*5, - 'X_h2': 245.67774602566578*5, - 'X_I': 206.42934561053158, - 'S_cat': 40.0, - 'S_an': 20.0, - } - -R1_ss_conds = { - 'S_su': 0.0145871088552909*1e3, - 'S_aa': 0.00643308564144693*1e3, - 'S_fa': 0.634823005990967*1e3, - 'S_va': 0.624510322247682*1e3, - 'S_bu': 1.03793927591996*1e3, - 'S_pro': 1.24676871525373*1e3, - 'S_ac': 2.00250371674824*1e3, - 'S_h2': 0.00850364943532684*1e3, - 'S_ch4': 0.0000422133982597226*1e3, - 'S_IC': 0.0951*C_mw*1e3, - 'S_IN': 0.0945*N_mw*1e3, - 'S_I': 0.027310256066728*1e3, - 'X_c': 0.146203507058736*1e3, - 'X_ch': 0.0286018513139117*1e3, - 'X_pr': 0.0467836694957302*1e3, - 'X_li': 0.0247209587890493*1e3, - 'X_su': 4.69052782535406*1e3, - 'X_aa': 1.22829926704024*1e3, - 'X_fa': 0.0147446263753011*1e3, - 'X_c4': 0.0149933579422897*1e3, - 'X_pro': 0.0145343147735253*1e3, - 'X_ac': 0.00098041337766024*1e3, - 'X_h2': 0.00110808891184369*1e3, - 'X_I': 0.0396205121367899*1e3 - } - -R2_ss_conds = { - 'S_su': 0.00106990968535691*1e3, - 'S_aa': 0.00125571416517827*1e3, - 'S_fa': 0.121097573221394*1e3, - 'S_va': 0.0132519103137696*1e3, - 'S_bu': 0.0172912281196732*1e3, - 'S_pro': 0.020032163173878*1e3, - 'S_ac': 0.00574002755366853*1e3, - 'S_h2': 3.76969944940856e-08*1e3, - 'S_ch4': 0.0499411746585487*1e3, - 'S_IC': 0.0951*C_mw*1e3, - 'S_IN': 0.0945*N_mw*1e3, - 'S_I': 0.105601391746794*1e3, - 'X_c': 0.0897520281015078*1e3, - 'X_ch': 0.00108163641708242*1e3, - 'X_pr': 0.00120204580901502*1e3, - 'X_li': 0.00150204523369107*1e3, - 'X_su': 0.195961987850137*1e3, - 'X_aa': 0.059723477130333*1e3, - 'X_fa': 0.0351858744892462*1e3, - 'X_c4': 0.0812315951844566*1e3, - 'X_pro': 0.0503466475437059*1e3, - 'X_ac': 1.1653549028287*1e3, - 'X_h2': 0.4352809013846*1e3, - 'X_I': 0.196117291164614*1e3 - } - -#%% Systems -if not _impact_item_loaded: load_lca_data() - -def create_systems(lifetime=30, discount_rate=0.1, electric_price=0.0913, - flowsheet_A=None, flowsheet_B=None, flowsheet_C=None, flowsheet_D=None, - inf_concs={}, R1_init_conds={}, R2_init_conds={}, R_init_conds={}, - which=None, selective=False): - PowerUtility.price = electric_price - which = which or ('A', 'B', 'C', 'D') - if isinstance(which, str): which = (which,) - - if selective: - R1_retain = fermenters - R2_retain = methanogens - else: R1_retain = R2_retain = biomass_IDs - pc.create_adm1_cmps() - adm1 = pc.ADM1() - adm1_phctrl = pc.ADM1() - dyn_params = adm1_phctrl.rate_function.params.copy() - adm1_phctrl.set_rate_function(rhos_adm1_ph_ctrl) - adm1_phctrl.rate_function._params = dyn_params - - inf_concs = inf_concs or default_inf_concs.copy() - brewery_ww = WasteStream('BreweryWW_A') - brewery_ww.set_flow_by_concentration(Q, concentrations=inf_concs, - units=('m3/d', 'kg/m3')) - R1_init_conds = R1_init_conds or R1_ss_conds - R2_init_conds = R2_init_conds or R2_ss_conds - R_init_conds = R_init_conds or default_R_init_conds - systems = [] - - if 'A' in which: - flowsheet_A = flowsheet_A or qs.Flowsheet('A') - qs.main_flowsheet.set_flowsheet(flowsheet_A) - flowsheet_A.stream.register(brewery_ww.ID, brewery_ww) - - effA = WasteStream('Effluent_A', T=T2) - bg1A = WasteStream('Biogas_1A', phase='g') - bg2A = WasteStream('Biogas_2A', phase='g') - - ISTA = IST(ID='ISTA') - GHA = GH(ID='GHA') - - R1A = AB('R1A', ins=brewery_ww, outs=(bg1A, ''), - V_liq=Vl1, V_gas=Vg1, T=T1, model=adm1_phctrl, - retain_cmps=R1_retain, bead_lifetime=bl, - exogenous_vars=(Temp1, pH1), - F_BM_default=1) - R2A = AB('R2A', ins=R1A-1, outs=(bg2A, effA), - V_liq=Vl2, V_gas=Vg2, T=T2, model=adm1_phctrl, - retain_cmps=R2_retain, bead_lifetime=bl, - exogenous_vars=(Temp2, pH2), - equipment=[ISTA, GHA], - F_BM_default=1) - - R1A.set_init_conc(**R1_init_conds) - R2A.set_init_conc(**R2_init_conds) - - sysA = System('sysA', path=(R1A, R2A)) - sysA.set_dynamic_tracker(R1A, R2A, bg1A, bg2A, effA) - add_strm_iitm(sysA) - add_TEA_LCA(sysA, discount_rate, lifetime) - systems.append(sysA) - - if 'B' in which: - flowsheet_B = flowsheet_B or qs.Flowsheet('B') - qs.main_flowsheet.set_flowsheet(flowsheet_B) - - infB = brewery_ww.copy('BreweryWW_B') - effB = WasteStream('Effluent_B', T=T2) - bg1B = WasteStream('Biogas_1B', phase='g') - bgh2B = WasteStream('Biogas_hsp_2B', phase='g') - bgm2B = WasteStream('Biogas_mem_2B', phase='g') - - ISTB = IST(ID='ISTB') - GHB = GH(ID='GHB') - - R1B = AB('R1B', ins=infB, outs=(bg1B, ''), - V_liq=Vl1, V_gas=Vg1, T=T1, model=adm1_phctrl, - retain_cmps=R1_retain, bead_lifetime=bl, - exogenous_vars=(Temp1, pH1), - F_BM_default=1) - R2B = AB('R2B', ins=R1B-1, outs=(bgh2B, ''), - V_liq=Vl2, V_gas=Vg2, T=T2, model=adm1_phctrl, - retain_cmps=R2_retain, bead_lifetime=bl, - exogenous_vars=(Temp2, pH2), - equipment=[ISTB, GHB], - F_BM_default=1) - DM2B = DM('DM2B', ins=R2B-1, outs=(bgm2B, effB), - F_BM_default=1) - - R1B.set_init_conc(**R1_init_conds) - R2B.set_init_conc(**R2_init_conds) - - sysB = System('sysB', path=(R1B, R2B, DM2B)) - sysB.set_dynamic_tracker(R1B, R2B, bg1B, bgh2B, bgm2B, effB) - add_strm_iitm(sysB) - add_TEA_LCA(sysB, discount_rate, lifetime) - systems.append(sysB) - - if 'C' in which: - flowsheet_C = flowsheet_C or qs.Flowsheet('C') - qs.main_flowsheet.set_flowsheet(flowsheet_C) - - infC = brewery_ww.copy('BreweryWW_C') - effC = WasteStream('Effluent_C', T=T2) - bgC = WasteStream('Biogas_C', phase='g') - - ISTC = IST(ID='ISTC') - GHC = GH(ID='GHC') - - RC = AB('RC', ins=infC, outs=(bgC, effC), - V_liq=Vl3, V_gas=Vg3, T=T3, model=adm1, - retain_cmps=biomass_IDs, bead_lifetime=bl, - equipment=[ISTC, GHC], - F_BM_default=1) - RC.set_init_conc(**R_init_conds) - - sysC = System('sysC', path=(RC,)) - sysC.set_dynamic_tracker(RC, bgC, effC) - add_strm_iitm(sysC) - add_TEA_LCA(sysC, discount_rate, lifetime) - systems.append(sysC) - - if 'D' in which: - flowsheet_D = flowsheet_D or qs.Flowsheet('D') - qs.main_flowsheet.set_flowsheet(flowsheet_D) - - infD = brewery_ww.copy('BreweryWW_D') - effD = WasteStream('Effluent_D', T=T2) - bghD = WasteStream('Biogas_hsp_D', phase='g') - bgmD = WasteStream('Biogas_mem_D', phase='g') - - ISTD = IST(ID='ISTD') - GHD = GH(ID='GHD') - - RD = AB('RD', ins=infD, outs=(bghD, ''), - V_liq=Vl3, V_gas=Vg3, T=T3, model=adm1, - retain_cmps=biomass_IDs, bead_lifetime=bl, - equipment=[ISTD, GHD], - F_BM_default=1) - DMD = DM('DMD', ins=RD-1, outs=(bgmD, effD), - F_BM_default=1) - RD.set_init_conc(**R_init_conds) - - sysD = System('sysD', path=(RD, DMD)) - sysD.set_dynamic_tracker(RD, bghD, bgmD, effD) - add_strm_iitm(sysD) - add_TEA_LCA(sysD, discount_rate, lifetime) - systems.append(sysD) - - return systems - -#%% -if __name__ == '__main__': - systems = create_systems() - for sys in systems: - sys.simulate( - t_span=(0,400), - state_reset_hook='reset_cache', - method='BDF' - ) - sys.diagram() \ No newline at end of file diff --git a/exposan/metab_mock/test_biofilm.py b/exposan/metab_mock/test_biofilm.py deleted file mode 100644 index e7883df0..00000000 --- a/exposan/metab_mock/test_biofilm.py +++ /dev/null @@ -1,598 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 10 15:53:18 2023 - -@author: joy_c -""" -import numpy as np, pandas as pd -import qsdsan.processes as pc -from qsdsan.utils import auom, ospath, time_printer -from scipy.integrate import solve_ivp -from exposan.metab_mock import results_path, biomass_IDs, create_systems - -Q = 5 # m3/d - -cmps = pc.create_adm1_cmps() -bm_idx = cmps.indices(biomass_IDs) -S_idx = cmps.indices([i for i in cmps.IDs if i.startswith('S_')]) -adm1 = pc.ADM1() -stoi_bk = adm1.stoichio_eval() -stoi_en = stoi_bk[:-3] # no liquid-gas transfer -Rho_bk = adm1.rate_function - -#%% functions & constants - -from qsdsan.processes._adm1 import ( - mass2mol_conversion, - T_correction_factor, - acid_base_rxn, - substr_inhibit, - Hill_inhibit, - non_compet_inhibit - ) -from scipy.optimize import brenth - -rhos = np.zeros(19) # 19 kinetic processes, no gas-liquid transfer -Cs = np.empty(19) -def Rho_en(state_arr, params=Rho_bk._params, T_op=298.15): - ks = params['rate_constants'] - Ks = params['half_sat_coeffs'] - cmps = params['components'] - pH_ULs = params['pH_ULs'] - pH_LLs = params['pH_LLs'] - KS_IN = params['KS_IN'] - KI_nh3 = params['KI_nh3'] - KIs_h2 = params['KIs_h2'] - Kab = params['Ka_base'] - Ka_dH = params['Ka_dH'] - T_base = params['T_base'] - - Cs[:8] = state_arr[12:20] - Cs[8:12] = state_arr[19:23] - Cs[12:] = state_arr[16:23] - - substrates = state_arr[:8] - - S_va, S_bu, S_h2, S_IN = state_arr[[3,4,7,10]] - unit_conversion = mass2mol_conversion(cmps) - cmps_in_M = state_arr[:27] * unit_conversion - weak_acids = cmps_in_M[[24, 25, 10, 9, 6, 5, 4, 3]] - - Kas = Kab * T_correction_factor(T_base, T_op, Ka_dH) - - rhos[:] = ks * Cs - Monod = substr_inhibit(substrates, Ks) - rhos[4:12] *= Monod - if S_va > 0: rhos[7] *= 1/(1+S_bu/S_va) - if S_bu > 0: rhos[8] *= 1/(1+S_va/S_bu) - - h = brenth(acid_base_rxn, 1e-14, 1.0, - args=(weak_acids, Kas), - xtol=1e-12, maxiter=100) - - nh3 = Kas[1] * weak_acids[2] / (Kas[1] + h) - - Iph = Hill_inhibit(h, pH_ULs, pH_LLs) - Iin = substr_inhibit(S_IN, KS_IN) - Ih2 = non_compet_inhibit(S_h2, KIs_h2) - Inh3 = non_compet_inhibit(nh3, KI_nh3) - rhos[4:12] *= Iph * Iin - rhos[6:10] *= Ih2 - rhos[10] *= Inh3 - return rhos - -n_cmps = len(cmps) -Diff = np.zeros(n_cmps) -Diff[:12] = [4.56e-6, 8.62e-6, 5.33e-6, 5.00e-6, 5.04e-6, 6.00e-6, - 6.48e-6, 4.02e-4, 1.36e-4, 1.71e-4, 1.52e-4, 6.00e-6] - -Diff[-3:-1] = 1.17e-4 -# f_diff = 0.75 # bead-to-water diffusivity fraction -Cs_in = np.array([ - 3.000e+03, 6.000e+02, 4.000e+02, 4.000e+02, 4.000e+02, 4.000e+02, - 4.000e+02, 5.000e-06, 5.000e-03, 4.804e+02, 1.401e+02, 2.000e+01, - 1.000e+02, 3.000e+02, 5.000e+02, 2.500e+02, 0.000e+00, 1.000e+00, - 1.000e+00, 1.000e+00, 1.000e+00, 1.000e+00, 1.000e+00, 2.500e+01, - 4.000e+01, 2.000e+01, 9.900e+05 - ])*1e-3 # mg/L to kg/m3 - -_R = 8.3145e-2 # Universal gas constant, [bar/M/K] -_k_p = 5.0e4 -_P_atm = 1.013 -Pa_to_bar = auom('Pa').conversion_factor('bar') -p_vapor = lambda T: cmps.H2O.Psat(T)*Pa_to_bar -gas_mass2mol_conversion = (cmps.i_mass / cmps.chem_MW)[cmps.indices(('S_h2', 'S_ch4', 'S_IC'))] -# K_tss = 25/2 # TSS half saturation coefficient - -def f_qgas(S_gas, T): - p_gas = S_gas * _R * T - P = sum(p_gas) + p_vapor(T) - q_gas = max(0, _k_p * (P - _P_atm)) - return q_gas - -def compile_ode(r_beads=5e-3, l_bl=1e-5, f_void=0.1, - n_dz=10, f_diff=0.75, K_tss=11): - dz = r_beads / n_dz - zs = np.linspace(dz, r_beads, n_dz) - dV = 4/3*np.pi*(zs)**3 - dV[1:] -= dV[:-1] - V_bead = (4/3*np.pi*r_beads**3) - D = Diff * f_diff # Diffusivity in beads - k = Diff/l_bl # mass transfer coeffient through liquid boundary layer - D_ov_dz2 = D[S_idx]/(dz**2) # (n_soluble,) - D_ov_dz = D[S_idx]/dz - _1_ov_z = 1/zs # (n_dz,) - - def dydt(t, y, HRT, detach=False): - Q, T = y[-2:] - V_liq = HRT * Q - V_bed = V_liq / f_void - V_beads = V_bed - V_liq - V_gas = V_bed * 0.1 - A_beads = 3 * V_beads / r_beads # m2, total bead surface area - - S_gas = y[-5:-2] - Cs_bk = y[-(n_cmps+5):-5] # bulk liquid concentrations - Cs_en = y[:n_dz*n_cmps].reshape((n_dz, n_cmps)) # each row is one control volume - - # Transformation - rhos_en = np.apply_along_axis(Rho_en, 1, Cs_en, T_op=T) - Rs_en = rhos_en @ stoi_en # n_dz * n_cmps - rhos_bk = Rho_bk(y[-(n_cmps+5):]) - Rs_bk = np.dot(stoi_bk.T, rhos_bk) # n_cmps+5 - q_gas = f_qgas(S_gas, T) - gas_transfer = - q_gas*S_gas/V_gas + rhos_bk[-3:] * V_liq/V_gas * gas_mass2mol_conversion - - # Detachment -- particulates - if detach: - tss = np.sum(Cs_en * (cmps.x*cmps.i_mass), axis=1) - x_net_growth = np.sum(Rs_en * cmps.x, axis=1)/np.sum(Cs_en * cmps.x, axis=1) # d^(-1), equivalent to k_de - u_de = 1/(1+np.exp(K_tss-tss)) * np.maximum(x_net_growth, 0) - de_en = np.diag(u_de) @ (Cs_en * cmps.x) - tot_de = np.sum(np.diag(dV) @ de_en, axis=0) / V_bead # detachment per unit volume of beads - else: - de_en = tot_de = 0 - - # Mass transfer (centered differences) -- MATLAB/Simulink - # S_en = Cs_en[:, S_idx] - # dSdz2 = np.zeros_like(Cs_en) - # dSdz2[0, S_idx] = (1+dz/(2*zs[0]))**2 * (S_en[1] - S_en[0])/dz - # dSdz2[1:-1, S_idx] = (np.diag((1+dz/(2*zs[1:-1]))**2) @ (S_en[2:]-S_en[1:-1]) \ - # - np.diag((1-dz/(2*zs[1:-1]))**2) @ (S_en[1:-1]-S_en[:-2]))/(dz**2) - # C_lf = Cs_en[-1].copy() - # S_en[-1] = C_lf[S_idx] = (D/dz*Cs_en[-2] + k*Cs_bk)[S_idx]/(D/dz+k)[S_idx] # Only solubles based on outer b.c. - # dSdz2[-1, S_idx] = ((1+dz/(2*zs[-1]))**2 * (S_en[-1]-S_en[-2]) \ - # -(1-dz/(2*zs[-1]))**2 * (S_en[-2]-S_en[-3]))/(dz**2) - - # dCdz = Cs_en * 0. - # d2Cdz2 = Cs_en * 0. - # dCdz[1:-1] = (Cs_en[2:] - Cs_en[:-2])/(2*dz) - # d2Cdz2[1:-1] = (Cs_en[2:] + Cs_en[:-2] - 2*Cs_en[1:-1])/(dz**2) - # d2Cdz2[0,:] = (Cs_en[1] - Cs_en[0])*(1+dz/(2*zs[0]))**2/dz # inner boundary condition, forward difference - # C_lf = Cs_en[-1].copy() - # C_lf[S_idx] = (D*Cs_en[-2] + k*Cs_bk*dz)[S_idx]/(D+k*dz)[S_idx] # Only solubles based on outer b.c. - # d2Cdz2[-1,:] = ((1+dz/(2*zs[-1]))**2*(C_lf-Cs_en[-2]) - (1-dz/(2*zs[-1]))**2*(Cs_en[-2]-Cs_en[-3]))/(dz**2) # backward difference - - #!!! Mass transfer (centered differences) -- MOL; solubles only - C_lf = Cs_en[-1] - J_lf = k*(Cs_bk - C_lf) - S_en = Cs_en[:, S_idx] - M_transfer = np.zeros_like(Cs_en) - M_transfer[1:-1, S_idx] = D_ov_dz2 * (S_en[2:] - 2*S_en[1:-1] + S_en[:-2])\ - + D_ov_dz * (np.diag(_1_ov_z[1:-1]) @ (S_en[2:] - S_en[:-2])) - M_transfer[0, S_idx] = 2 * D_ov_dz2 * (S_en[1] - S_en[0]) - M_transfer[-1, S_idx] = 2 * D_ov_dz2 * (S_en[-2] - S_en[-1])\ - + 2 * (1/dz + _1_ov_z[-1]) * J_lf[S_idx] - - # Mass balance - # dCdt_en = D*dSdz2 + Rs_en - de_en - # dCdt_en = D*(d2Cdz2 + np.diag(2/zs) @ dCdz) + Rs_en - de_en - # dCdt_en = dJ_dz + Rs_en - de_en - dCdt_en = M_transfer + Rs_en - de_en - - # dCdt_bk = (Cs_in-Cs_bk)/HRT - A_beads/V_liq*k*(Cs_bk-C_lf) + Rs_bk + V_beads*tot_de/V_liq - dCdt_bk = (Cs_in-Cs_bk)/HRT - A_beads/V_liq*J_lf + Rs_bk + V_beads*tot_de/V_liq - - dy = np.zeros_like(y) - dy[:n_dz*n_cmps] = dCdt_en.flatten() - dy[-(n_cmps+5):-5] = dCdt_bk - dy[-5:-2] = gas_transfer - - return dy - - return dydt - -#%% -# HRT = 1 -C0 = np.array([ - 1.204e-02, 5.323e-03, 9.959e-02, 1.084e-02, 1.411e-02, 1.664e-02, - 4.592e-02, 2.409e-07, 7.665e-02, 5.693e-01, 1.830e-01, 3.212e-02, - 2.424e-01, 2.948e-02, 4.766e-02, 2.603e-02, 4.708e+00, 1.239e+00, - 4.838e-01, 1.423e+00, 8.978e-01, 2.959e+00, 1.467e+00, 4.924e-02, - 4.000e-02, 2.000e-02, 9.900e+02 - ]) - -y0_bulk = np.array([ - 1.204e-02, 5.323e-03, 9.959e-02, 1.084e-02, 1.411e-02, 1.664e-02, - 4.592e-02, 2.409e-07, 7.665e-02, 5.693e-01, 1.830e-01, 3.212e-02, - 2.424e-01, 2.948e-02, 4.766e-02, 2.603e-02, 0, 0, - 0, 0, 0, 0, 0, 0, # no biomass in bulk liquid initially - 4.000e-02, 2.000e-02, 9.900e+02, 3.922e-07, 2.228e-02, 1.732e-02, - Q, 298.15 - ]) - -y0_from_bulk = lambda n_dz: np.append(np.tile(C0, n_dz), y0_bulk) - -#%% steady-state spatial profile -def y0_even(n_dz, TSS0): - C0_even = C0.copy() - # C0_even = Cs_in.copy() - C0_even[bm_idx] = TSS0/0.777/7 - return np.append(np.tile(C0_even, n_dz), y0_bulk) - -def spatial_profiling(HRTs=[1, 0.5, 10/24, 8/24], TSS0=5, n_dz=10, - detach=False, save_to='', **ode_kwargs): - y0 = y0_even(n_dz, TSS0) - # y0 = y0_from_bulk(n_dz) - dfs = {} - dydt = compile_ode(n_dz=n_dz, **ode_kwargs) - for tau in HRTs: - print(f'HRT = {tau:.2f} d') - sol = solve_ivp(dydt, t_span=(0, 400), y0=y0, method='BDF', args=(tau, detach)) - y_ss = sol.y.T[-1] - C_ss = y_ss[:n_cmps*(n_dz+1)].reshape((n_dz+1, n_cmps)) - Xbio_ss = C_ss[:, bm_idx] - df_c = pd.DataFrame(C_ss, columns=cmps.IDs, index=[*range(n_dz), 'bulk']) - df_c['biomass_COD'] = np.sum(Xbio_ss, axis=1) - dfs[f'{tau:.2f}'] = df_c - - file_name = save_to or f'TSS0_{TSS0}_{detach}.xlsx' - path = ospath.join(results_path, file_name) - with pd.ExcelWriter(path) as writer: - for k, df in dfs.items(): - df.to_excel(writer, sheet_name=k) - return dfs - -#%% suspended vs. attached growth - -@time_printer -def suspended_vs_attached(HRT, dydt, n_dz, r_beads, sys, f_void=0.1, - y0=None, detach=False): - print('\n') - print('='*12) - print(f'HRT: {HRT:.3f} d') - u, = sys.units - u.V_liq = HRT * 5 - u.V_gas = HRT * 0.5 - inf, = sys.feeds - bg, eff = sys.products - sys.simulate(state_reset_hook='reset_cache', t_span=(0, 400), method='BDF') - sus_rcod = 1-eff.COD/inf.COD - i = 0 - while sus_rcod < 0.7 and i <= 5: - i += 1 - print('add VSS...') - u._concs[bm_idx] *= 2 - sys.simulate(state_reset_hook='reset_cache', t_span=(0, 400), method='BDF') - sus_rcod = 1-eff.COD/inf.COD - sus_bm = sum(u._state[bm_idx] * cmps.i_mass[bm_idx]) - if y0 is None: - y0_susp = np.append(u._state, 298.15) - y0_susp[bm_idx] = 0. - y0 = np.append(np.tile(u._state[:n_cmps], n_dz), y0_susp) - try: y = converge(dydt, y0, args=(HRT, detach)) - except: - try: - y0 = y0_even(n_dz, TSS0=min(2.5/HRT, 22)) - y = converge(dydt, y0, args=(HRT, detach)) - except: - return sus_rcod, np.nan, sus_bm, np.nan, None - bk_Cs = y[-(n_cmps+5):-5].copy() - att_COD = sum(bk_Cs * cmps.i_COD * (1-cmps.g)) * 1e3 - att_rcod = 1-att_COD/inf.COD - bk_bm, en_bm = biomass_CODs(y, n_dz, r_beads) - att_bm = ((1-f_void)*en_bm + f_void*bk_bm) * 0.777 - print(bk_bm*0.777, en_bm*0.777) - return [sus_rcod, att_rcod, sus_bm, att_bm, y] - -def HRT_suspended_vs_encap(HRTs=[1, 0.5, 10/24, 8/24, 4/24, 2/24, 1/24], - frac_retain=1.0, r_beads=1e-5, l_bl=1e-6, n_dz=8, - **ode_kwargs): - record = [] - dydt = compile_ode(r_beads=r_beads, l_bl=l_bl, n_dz=n_dz, **ode_kwargs) - - sys, = create_systems(which='C') - u, = sys.units - u._f_retain = (u._f_retain > 0) * frac_retain - u.encapsulate_concentration = 1e5 - y = y0_even(n_dz, TSS0=2.5) - # y = None - for tau in HRTs: - out = suspended_vs_attached(tau, dydt, n_dz, r_beads, sys, y, - detach=True) - y = out.pop(-1) - record.append(out) - - df = pd.DataFrame(record, index=HRTs, - columns=pd.MultiIndex.from_product( - [['rCOD', 'biomass [g TSS/L]'], ['Suspended', 'Encapsulated']] - )) - df.index.name = 'HRT [d]' - df.to_excel(ospath.join(results_path, 'suspended_vs_encap_retain.xlsx')) - return df - -#%% encapsulated biomass, HRT -def biomass_CODs(y, n_dz, r_beads): - en_bm = np.sum(y[:n_dz*n_cmps].reshape((n_dz, n_cmps))[:,bm_idx], axis=1) - bk_bm = np.sum(y[n_dz*n_cmps: ((n_dz+1)*n_cmps)][bm_idx]) - dz = r_beads / n_dz - zs = np.linspace(dz, r_beads, n_dz) - dV = 4/3*np.pi*(zs)**3 - V_bead = dV[-1] - dV[1:] -= dV[:-1] - C_en_avg = np.dot(en_bm, dV)/V_bead - return bk_bm, C_en_avg - -def SRT(y, n_dz, HRT, r_beads, f_void): - '''SRT given there's no retention of biomass in bulk liquid''' - bk_bm, C_en_avg = biomass_CODs(y, n_dz, r_beads) - return (1 + (1-f_void)/f_void*C_en_avg/bk_bm) * HRT - -def converge(dydt, y, args, threshold=1e-4, once=False): - dy_max = 1. - while dy_max > threshold: - sol = solve_ivp(dydt, t_span=(0, 400), y0=y, method='BDF', args=args) - y = sol.y.T[-1] - rcod = 1 - sum(y[-(n_cmps+5):-5] * cmps.i_COD * (1-cmps.g))/6.801 - dy = dydt(0, y, *args) - dy_max = np.max(np.abs(dy)) - print(f'rCOD {rcod:.3f}, dy_max {dy_max:.2e}') - if once: break - return y - -@time_printer -def encap(TSS_init, HRT, dydt, n_dz, f_void=0.1, r_beads=5e-3, - detach=True, y0=None, threshold=1e-4): - msg = f'r_beads {r_beads*1e3:.2f}mm, HRT {HRT:.2f}d' - print(f'\n{msg}\n{"="*len(msg)}') - if y0 is None: y = y0_even(n_dz, TSS_init) - else: y = y0 - try: - y = converge(dydt, y, args=(HRT, detach), threshold=threshold) - except: - try: - print('Try bulk y0...') - # y = y0_even(n_dz, TSS_init*2) - y = y0_from_bulk(n_dz) - y = converge(dydt, y, args=(HRT, detach), threshold=threshold) - except: return [np.nan]*4 - bk_Cs = y[-(n_cmps+5):-5] - en_Cs = y[-(2*n_cmps+5):-(n_cmps+5)] - att_COD = sum(bk_Cs * cmps.i_COD * (1-cmps.g)) * 1e3 - att_rcod = 1-att_COD/6801 - att_bm = sum(en_Cs[bm_idx] * cmps.i_mass[bm_idx]) - srt = SRT(y, n_dz, HRT, r_beads, f_void) - return att_rcod, att_bm, srt, y - -def HRT_init_TSS(TSS=[1, 2, 5, 10, 30], - HRTs=[1, 0.5, 10/24, 8/24, 4/24, 2/24, 1/24], - detach=False, **ode_kwargs): - rcod = [] - ss_bm = [] - dydt = compile_ode(detach=detach, **ode_kwargs) - r_beads = ode_kwargs.pop('r_beads', 5e-3) - f_void = ode_kwargs.pop('f_void', 0.1) - n_dz = ode_kwargs.pop('n_dz', 10) - for tss in TSS: - l_rcod = [] - l_bm = [] - for tau in HRTs: - out = encap(tss, tau, dydt, n_dz, f_void=f_void, r_beads=r_beads, - detach=detach) - l_rcod.append(out[0]) - l_bm.append(out[1]) - rcod.append(l_rcod) - ss_bm.append(l_bm) - - df_rcod = pd.DataFrame(rcod).T - df_bm = pd.DataFrame(ss_bm).T - - for i in (df_rcod, df_bm): - i.columns = pd.MultiIndex.from_product([['initial TSS [g/L]'], TSS]) - i.index = HRTs - i.index.name = 'HRT [d]' - - with pd.ExcelWriter(ospath.join(results_path, 'HRT_vs_encapTSS_detach.xlsx')) as writer: - df_rcod.to_excel(writer, sheet_name='rCOD') - df_bm.to_excel(writer, sheet_name='Biomass TSS') - -#%% bead size - -def bead_size_HRT(HRTs=np.array([1, 0.5, 10/24, 8/24, 4/24, 2/24, 1/24]), - bead_size=np.array([5e-3, 1e-3, 5e-4, 1e-4, 1e-5]), - voidage=[0.39, 0.4, 0.5, 0.6, 0.7], - # voidage=[0.39]*5, - detach=False, **ode_kwargs): - rcod = [] - ss_bm = [] - srt = [] - for r_beads, f_void in zip(bead_size, voidage): - l_bl = min(1e-5, r_beads/10) - if r_beads > 1e-3: n_dz = 10 - elif r_beads < 1e-4: n_dz = 5 - else: n_dz = 8 - dydt = compile_ode(r_beads=r_beads, l_bl=l_bl, f_void=f_void, n_dz=n_dz, **ode_kwargs) - l_rcod = [] - l_bm = [] - l_srt = [] - y0 = y0_from_bulk(n_dz) - for tau in HRTs: - tss = min(2.5/tau, 22) - z1, z2, z3, y0 = encap(tss, tau, dydt, n_dz, f_void, - r_beads, detach, y0, 5e-4) - l_rcod.append(z1) - l_bm.append(z2) - l_srt.append(z3) - rcod.append(l_rcod) - ss_bm.append(l_bm) - srt.append(l_srt) - - # sys, = create_systems(which='C') - # u, = sys.units - # inf, = sys.feeds - # bg, eff = sys.products - # u.encapsulate_concentration = 1e5 - # l_rcod = [] - # for tau, tau_s in zip(HRTs, srt[-1]): - # u._f_retain[bm_idx] = 1-tau/tau_s - # u.V_liq = tau*5 - # u.V_gas = tau*0.5 - # sys.simulate(state_reset_hook='reset_cache', - # t_span=(0, 400), method='BDF') - # l_rcod.append(1-eff.COD/inf.COD) - - df_rcod = pd.DataFrame(rcod).T - df_bm = pd.DataFrame(ss_bm).T - df_srt = pd.DataFrame(srt).T - - for i in (df_rcod, df_bm, df_srt): - # i.columns = pd.MultiIndex.from_product([['bead size [m]'], bead_size]) - i.columns = pd.MultiIndex.from_tuples(zip(np.round(bead_size*1e3, 2), voidage)) - i.columns.names = ['bead size [mm]', 'voidage'] - i.index = np.int32(HRTs*24) - i.index.name = 'HRT [h]' - # df_rcod[('Suspended','')] = l_rcod - with pd.ExcelWriter(ospath.join(results_path, - 'HRT_vs_rbeads_incre.xlsx')) as writer: - df_rcod.to_excel(writer, sheet_name='rCOD') - df_bm.to_excel(writer, sheet_name='encapsulated TSS') - df_srt.to_excel(writer, sheet_name='SRT') - - return df_rcod, df_bm, df_srt - -#%% -def voidage_vs_hrt(voidage = np.linspace(0.39, 0.8, 10), - HRTs=np.array([1, 0.5, 10/24, 8/24, 4/24, 2/24, 1/24]), - r_beads=1e-5, l_bl=1e-6, n_dz=5, f_diff=0.75, K_tss=11): - sys, = create_systems(which='C') - u, = sys.units - inf, = sys.feeds - bg, eff = sys.products - u.encapsulate_concentration = 1e5 - # df_rcod = pd.DataFrame(index=np.int32(HRTs*24), - df_rcod = pd.DataFrame(index=np.round(HRTs*24, 1), - columns=pd.MultiIndex.from_product( - [np.round(voidage,3), ['encap','suspend']] - )) - df_rcod.index.name = 'HRT [h]' - df_rcod.columns.names = ['voidage', 'biomass'] - df_bm = df_rcod.copy() - df_srt = df_rcod.copy() - # df_srt = pd.DataFrame(index=np.int32(HRTs*24), columns=np.round(voidage,3)) - # df_srt.index.name = df_rcod.index.name - # df_srt.columns.name = 'voidage' - for f_void in voidage: - dydt = compile_ode(r_beads=r_beads, l_bl=l_bl, f_void=f_void, n_dz=n_dz, - f_diff=f_diff, K_tss=K_tss) - y0 = y0_from_bulk(n_dz) - y = y0.copy() - for tau in HRTs: - # itau = int(tau*24) - itau = round(tau*24, 1) - ivoid = round(f_void,3) - msg = f'voidage {ivoid}, HRT {itau}h' - print(f'\n{msg}\n{"="*len(msg)}') - if f_void < 0.4: y = y0.copy() - y = converge(dydt, y, (tau, True), - threshold=5e-4, once=False) - bk_Cs = y[-(n_cmps+5):-5].copy() - att_COD = sum(bk_Cs * cmps.i_COD * (1-cmps.g)) * 1e3 - df_rcod.loc[itau, (ivoid, 'encap')] = 1-att_COD/inf.COD - bk_bm, en_bm = biomass_CODs(y, n_dz, r_beads) - df_bm.loc[itau, (ivoid, 'encap')] = ((1-f_void)*en_bm + f_void*bk_bm) * 0.777 - df_srt.loc[itau, (ivoid, 'encap')] = srt = SRT(y, n_dz, tau, r_beads, f_void) - - # u._f_retain[bm_idx] = 1-tau/srt - u._f_retain[:] = (1-tau/srt) * cmps.x - u.V_liq = tau*5 - u.V_gas = tau*0.5 - sys.simulate(state_reset_hook='reset_cache', - t_span=(0, 400), method='BDF') - df_rcod.loc[itau, (ivoid, 'suspend')] = 1-eff.COD/inf.COD - df_bm.loc[itau, (ivoid, 'suspend')] = sp_bm = sum(u._state[bm_idx]) * 0.777 - df_srt.loc[itau, (ivoid, 'suspend')] = tau * sp_bm/(sum(eff.state[bm_idx])/1e3*0.777) - - with pd.ExcelWriter(ospath.join(results_path, 'void_vs_HRT_opty0.xlsx')) as writer: - df_rcod.to_excel(writer, sheet_name='rCOD') - df_bm.to_excel(writer, sheet_name='biomass TSS') - df_srt.to_excel(writer, sheet_name='SRT') - - return df_rcod, df_bm, df_srt - -#%% -@time_printer -def run(tau=0.5, TSS0=5, r_beads=5e-3, l_bl=1e-5, f_void=0.1, n_dz=10, - f_diff=0.75, K_tss=11, detach=True, y0=None): - if y0 is None: - # y0 = y0_even(n_dz, TSS0) - y0 = y0_from_bulk(n_dz) - dydt = compile_ode(r_beads, l_bl, f_void, n_dz, f_diff, K_tss) - # print(f'HRT = {tau:.2f} d') - y_ss = converge(dydt, y0, args=(tau, detach), threshold=5e-4) - C_ss = y_ss[:n_cmps*(n_dz+1)].reshape((n_dz+1, n_cmps)) - Xbio_ss = C_ss[:, bm_idx] - df_c = pd.DataFrame(C_ss, columns=cmps.IDs, index=[*range(1, n_dz+1), 'bulk']) - df_c['biomass_COD'] = np.sum(Xbio_ss, axis=1) - bk = C_ss[-1].copy() - eff_cod = np.sum(bk * cmps.i_COD * (1-cmps.g)) - rcod = 1-eff_cod/6.801 - srt = SRT(y_ss, n_dz=n_dz, HRT=tau, r_beads=r_beads, f_void=f_void) - return df_c, y_ss, rcod, srt - - -if __name__ == '__main__': - # dfs = bead_size_HRT(detach=True) - # dfs = bead_size_HRT(bead_size=[1e-4, 1e-5], detach=True, - # n_dz=8, f_void=0.95) - # de_rcod, de_bm = bead_size_HRT(detach=True, n_dz=20) - - # for detach in (True,): - # print(f'detach: {detach}') - # for r_beads in (5e-3, 1e-3, 5e-4, 1e-4, 1e-5): - # print(f'r_beads: {r_beads}\n{"="*15}') - # dfs = spatial_profiling( - # HRTs=[1, 0.5, 1/3], TSS0=6, - # r_beads=r_beads, detach=detach, - # save_to=f'spatial_r{r_beads}_{detach}_new.xlsx' - # ) - # df, y, rcod, srt = run(tau=1, TSS0=2.5, r_beads=1e-5, l_bl=1e-6, n_dz=5, - # f_void=0.5, detach=True) - # df = run(tau=1/6, TSS0=6, detach=True) - - # df = HRT_suspended_vs_encap() - rcod, tss, srt = voidage_vs_hrt() - # dfs = {} - # rcods = [] - # srts = [] - # for f_void in (0.6, 0.8): - # l_rcod = [] - # l_srt = [] - # for n_dz in (5, 10, 20): - # print(f'\nvoidage {f_void}, n = {n_dz}') - # df, y, rcod, srt = run(r_beads=1e-5, l_bl=1e-6, - # f_void=f_void, n_dz=n_dz) - # dfs[f'void {f_void} (n={n_dz})'] = df - # l_rcod.append(rcod) - # l_srt.append(srt) - # rcods.append(l_rcod) - # srts.append(l_srt) - # summary = rcods + srts - # summary = pd.DataFrame(summary).T - # summary.index = [5,10,20] - # summary.index.name = 'n_dz' - # summary.columns = pd.MultiIndex.from_product([['rCOD', 'SRT [d]'], [0.6, 0.8]]) - # summary.columns.names = ['','voidage'] - # with pd.ExcelWriter(ospath.join(results_path, 'ndz_vs_spatial.xlsx')) as writer: - # summary.to_excel(writer, sheet_name='summary') - # for k,v in dfs.items(): - # v.to_excel(writer, sheet_name=k) \ No newline at end of file diff --git a/exposan/metab_mock/units.py b/exposan/metab_mock/units.py deleted file mode 100644 index 965f0f91..00000000 --- a/exposan/metab_mock/units.py +++ /dev/null @@ -1,303 +0,0 @@ -# -*- coding: utf-8 -*- -''' -EXPOsan: Exposition of sanitation and resource recovery systems - -This module is developed by: - - Joy Zhang - -This module is under the University of Illinois/NCSA Open Source License. -Please refer to https://github.com/QSD-Group/EXPOsan/blob/main/LICENSE.txt -for license details. -''' -from biosteam import Stream -from qsdsan import Construction -from qsdsan.sanunits import AnaerobicCSTR, Pump, HXutility -from qsdsan.utils import auom -from exposan.metab.equipment import Beads -from exposan.metab.utils import ( - pipe_design, pipe_friction_head, hdpe_price, - heat_transfer_U, - stainless_steel_wall_thickness as wt_ssteel, - add_prefix, - _construct_water_pump, - ) - -from math import pi -from collections import defaultdict - -__all__ = ('METAB_AnCSTR',) - - -#%% METAB_AnCSTR -_fts2mhr = auom('ft/s').conversion_factor('m/hr') -_cmph_2_gpm = auom('m3/hr').conversion_factor('gpm') - -class METAB_AnCSTR(AnaerobicCSTR): - - auxiliary_unit_names = ('heat_exchanger', ) - - def __init__(self, ID='', lifetime=30, bead_lifetime=10, - encapsulate_concentration=25, - reactor_height_to_diameter=1.5, - mixing_power=10, # hp/1000 gallons - wall_concrete_unit_cost=1081.73, # $850/m3 in 2014 USD, converted to 2021 USD with concrete PPI - slab_concrete_unit_cost=582.48, # $458/m3 in 2014 USD - stainless_steel_unit_cost=1.8, # https://www.alibaba.com/product-detail/brushed-stainless-steel-plate-304l-stainless_1600391656401.html?spm=a2700.details.0.0.230e67e6IKwwFd - rockwool_unit_cost=0.59, # https://www.alibaba.com/product-detail/mineral-wool-insulation-price-mineral-wool_60101640303.html?spm=a2700.7724857.0.0.262334d1rZXb48 - carbon_steel_unit_cost=0.5, # https://www.alibaba.com/product-detail/ASTM-A106-Ss400-Q235-Standard-Ms_1600406694387.html?s=p - **kwargs): - equip = kwargs.pop('equipment', []) - equip.append(Beads(ID=f'{ID}_beads', lifetime=bead_lifetime)) - super().__init__(ID, lifetime=lifetime, equipment=equip, **kwargs) - self.mixing_power=mixing_power - self.encapsulate_concentration = encapsulate_concentration - self.reactor_height_to_diameter = reactor_height_to_diameter - self.wall_concrete_unit_cost = wall_concrete_unit_cost - self.slab_concrete_unit_cost = slab_concrete_unit_cost - self.stainless_steel_unit_cost = stainless_steel_unit_cost - self.rockwool_unit_cost = rockwool_unit_cost - self.carbon_steel_unit_cost = carbon_steel_unit_cost - hx_in = Stream(f'{ID}_hx_in') - hx_out = Stream(f'{ID}_hx_out') - self.heat_exchanger = HXutility(ID=f'{ID}_hx', ins=hx_in, outs=hx_out) - for i in ('Wall concrete', 'Slab concrete', 'Stainless steel', - 'Rockwool', 'Carbon steel', 'HDPE pipes'): - name = i.lower().replace(' ', '_') - self.construction.append( - Construction(ID=name, linked_unit=self, item=name) - ) - for aux in self.auxiliary_units: - self.construction += aux.construction - for equip in self.equipment: - self.construction += equip.construction - - @property - def bead_lifetime(self): - for equip in self.equipment: - if isinstance(equip, Beads): return equip.lifetime - - @bead_lifetime.setter - def bead_lifetime(self, lt): - for equip in self.equipment: - if isinstance(equip, Beads): - equip.update_lifetime(lt) - - def _setup(self): - hasfield = hasattr - setfield = setattr - for i, ws in enumerate(self.ins): - field = f'Pump_ins{i}' - if not hasfield(self, field): - pump = Pump(ws.ID+'_Pump', ins=Stream(f'{ws.ID}_proxy')) - setfield(self, field, pump) - self.construction += [ - Construction(ID='22kW', linked_unit=pump, item='pump_22kW'), - Construction(ID='40W', linked_unit=pump, item='pump_40W') - ] - self.auxiliary_unit_names = tuple({*self.auxiliary_unit_names, field}) - super()._setup() - - _steel_insulate_thickness = 50 # mm - _concrete_cover_thickness = 100 # mm - _concrete_wall_thickness = 150 # mm - _concrete_base_thickness = 160 # mm - _cncr_insulate_thickness = 25 # mm - _facing_thickness = 3 # mm - - _density = { - 'Aluminum': 2710, # 2,640 - 2,810 kg/m3 - 'Stainless steel': 7930, # kg/m3, 18/8 Chromium - 'Rockwool': 100, # kg/m3 - 'Carbon steel': 7840 - } - - T_air = 273.15 + 20 - T_earth = 273.15 + 20 - - _l_min_velocity = 3*_fts2mhr - _g_min_velocity = 10*_fts2mhr - - _units = { - 'Volume': 'm3', - 'Height': 'm', - 'Outer diameter': 'm', - 'Area': 'm2', - 'Wall concrete': 'm3', - 'Slab concrete': 'm3', - 'Stainless steel': 'kg', - 'Rockwool': 'kg', - 'Carbon steel': 'kg', - 'HDPE pipes': 'kg', - 'Bead volume': 'm3', - 'Agitator': 'hp' - } - - def _design(self): - D = self.design_results - den = self._density - V = D['Volume'] = self.V_liq + self.V_gas - dia = (V*4/self.reactor_height_to_diameter/pi) ** (1/3) - h = D['Height'] = dia * self.reactor_height_to_diameter - tface = self._facing_thickness/1e3 - if V >= 5: - twall = self._concrete_wall_thickness/1e3 - tcover = self._concrete_cover_thickness/1e3 - tbase = self._concrete_base_thickness/1e3 - tinsl = self._cncr_insulate_thickness/1e3 - D['Outer diameter'] = OD = dia + twall * 2 - S_wall = pi*OD*h - S_base = D['Area'] = pi*(OD/2)**2 - D['Wall concrete'] = S_wall * twall - D['Slab concrete'] = S_base*(tcover + tbase) - D['Stainless steel'] = 0 - D['Rockwool'] = S_wall * tinsl * den['Rockwool'] - D['Carbon steel'] = S_wall * tface * den['Carbon steel'] - Uwall, Ucover, Ubase = heat_transfer_U(twall, tinsl, tface, tbase, concrete=True) - else: - P_liq = self.external_P * 1e5 + self._mixed.rho * 9.80665 * h * self.V_liq/V - twall = tbase = wt_ssteel(P_liq, dia, h) - tcover = twall/2 - tinsl = self._steel_insulate_thickness/1e3 - D['Outer diameter'] = OD = dia + twall*2 - S_wall = pi*OD*h - S_base = D['Area'] = pi*(OD/2)**2 - D['Wall concrete'] = 0 - D['Slab concrete'] = 0 - D['Stainless steel'] = (S_wall * twall + S_base*(tcover + tbase)) * den['Stainless steel'] - D['Rockwool'] = (S_base+S_wall) * tinsl * den['Rockwool'] - D['Carbon steel'] = (S_base+S_wall) * tface * den['Carbon steel'] - Uwall, Ucover, Ubase = heat_transfer_U(twall, tinsl, tface, tbase, concrete=False) - - # Calculate needed heating - T = self.T - hx = self.heat_exchanger - mixed = self._mixed - hx_ins0, hx_outs0 = hx.ins[0], hx.outs[0] - mixed.mix_from(self.ins) - hx_ins0.copy_flow(mixed) - hx_outs0.copy_flow(mixed) - hx_ins0.T = mixed.T - hx_outs0.T = T - hx_ins0.P = hx_outs0.P = mixed.T - - # Heat loss - wall_loss = Uwall * S_wall * (T-self.T_air) # [W] - base_loss = Ubase * S_base * (T-self.T_earth) # [W] - cover_loss = Ucover * S_base * (T-self.T_air) # [W] - duty = (wall_loss+base_loss+cover_loss)*60*60/1e3 # kJ/hr - hx.H = hx_ins0.H + duty # stream heating and heat loss - hx.simulate_as_auxiliary_exchanger(ins=hx.ins, outs=hx.outs) - if T <= mixed.T: hx.heat_utilities[0].cost = 0 - - # Piping - L_inlets = OD * 1.25 - L_outlets = h + OD*0.25 - L_gas = h + OD - pipe_IDs = [] - HDPE_pipes = [] - for ws in self.ins: - _inch, _kg_per_m = pipe_design(ws.F_vol, self._l_min_velocity) - pipe_IDs.append(_inch) - HDPE_pipes.append(_kg_per_m*L_inlets) - for ws in self.outs: - if ws.phase == 'g': - D['Stainless steel'] += pipe_design(ws.F_vol, self._g_min_velocity, True)[1] * L_gas - else: - _inch, _kg_per_m = pipe_design(ws.F_vol, self._l_min_velocity) - pipe_IDs.append(_inch) - HDPE_pipes.append(_kg_per_m*L_outlets) - D['HDPE pipes'] = sum(HDPE_pipes) - self._hdpe_ids, self._hdpe_kgs = pipe_IDs, HDPE_pipes - - # Vertical turbine agitator - D['Agitator'] = hp = self.mixing_power * auom('m3').convert(self.V_liq, 'gallon')/1000 - self.add_power_utility(auom('hp').convert(hp, 'kW')) - - # Pumps - flowsheet_ID = self.system.flowsheet.ID - getfield = getattr - creg = Construction.registry - for i, ws, in enumerate(self.ins): - ID = pipe_IDs[i] - hf = pipe_friction_head(ws.F_vol*_cmph_2_gpm, L_inlets, ID) # friction head loss - TDH = hf + h # in m, assume suction head = 0, discharge head = reactor height - field = f'Pump_ins{i}' - pump = getfield(self, field) - pump.ins[0].copy_flow(ws) - pump.dP_design = TDH * 9804.14 # in Pa - pump.simulate() - if self.include_construction: - q22, q40 = _construct_water_pump(pump) - p22 = getfield(creg, f'{flowsheet_ID}_{pump.ID}_22kW') - p40 = getfield(creg, f'{flowsheet_ID}_{pump.ID}_40W') - p22.quantity = q22 - p40.quantity = q40 - - if self.include_construction: - for i in ('Wall concrete', 'Slab concrete', 'Stainless steel', - 'Rockwool', 'Carbon steel', 'HDPE pipes'): - name = i.lower().replace(' ', '_') - const = getfield(creg, f'{flowsheet_ID}_{self.ID}_{name}') - const.quantity = D[i] - - # Beads - cmps = mixed.components - conc = self._state[:len(cmps)] - retained_concentration = sum(conc * cmps.i_mass * (self._f_retain > 0)) # kg mass/m3 - V_beads = D['Bead volume'] = self.V_liq * retained_concentration / self.encapsulate_concentration - if V_beads >= self.V_liq: - raise RuntimeError('retained biomass concentration > design encapsualation concentration') - self.add_equipment_design() - - _NG_price = 0.85*auom('kJ').conversion_factor('therm') # [USD/kJ] 5.47 2021USD/Mcf vs. 4.19 2016USD/Mcf - - def _cost(self): - bg = self.outs[0] - cmps = bg.components - KJ_per_kg = cmps.i_mass/cmps.chem_MW*cmps.LHV - self.add_OPEX['NG_offset'] = -sum(bg.mass*KJ_per_kg)*self._NG_price # kJ/hr * USD/kJ = USD/hr - D = self.design_results - C = self.baseline_purchase_costs - C['Wall concrete'] = D['Wall concrete']*self.wall_concrete_unit_cost - C['Slab concrete'] = D['Slab concrete']*self.slab_concrete_unit_cost - C['Stainless steel'] = D['Stainless steel']*self.stainless_steel_unit_cost - C['Rockwool'] = D['Rockwool']*self.rockwool_unit_cost - C['Carbon steel'] = D['Carbon steel']*self.carbon_steel_unit_cost - C['HDPE pipes'] = sum(hdpe_price(inch)*kg for inch, kg in zip(self._hdpe_ids, self._hdpe_kgs)) - hp = D['Agitator'] - D['Stainless steel'] += 29.792*hp**0.4785 # Empirical function based on data on https://www.agitadoresfluidmix.com/wp-content/uploads/2023/01/4.Data%20sheet%20VTA%20ENG.pdf - if hp < 2: C['Agitator'] = 2610*hp**0.34 * 708/567 # Seider et al., 2017, pp481; adjusted to 2021$ using CEPCI - else: C['Agitator'] = 3730*hp**0.54 * 708/567 - self.add_equipment_cost() - - def add_equipment_design(self): - unit_design = self.design_results - unit_units = self._units - isa = isinstance - get = getattr - if isa(self.equipment_lifetime, int): - lt = self.equipment_lifetime - self.equipment_lifetime = defaultdict(lambda: lt) - F_BM, F_D, F_P, F_M, lifetime = \ - self.F_BM, self.F_D, self.F_P, self.F_M, self.equipment_lifetime - - for equip in self.equipment: - equip_ID = equip.ID - prefix = f'{equip.__class__.__name__} {equip_ID}' - equip_design = equip._design_results = equip._design() - equip_design = {} if not equip_design else equip_design - unit_design.update(add_prefix(equip_design, prefix)) - - equip_units = {} if not equip.units else equip.units - unit_units.update(add_prefix(equip_units, prefix)) - for unit_attr, equip_attr in zip( - (F_BM, F_D, F_P, F_M, lifetime), - ('F_BM', 'F_D', 'F_P', 'F_M', 'lifetime'), - ): - equip_attr = get(equip, equip_attr) - if isa(equip_attr, dict): - unit_attr.update(add_prefix(equip_attr, prefix)) - else: - unit_attr[equip_ID] = equip_attr -