Skip to content

Commit

Permalink
MCAS w/generic crystallizer (watertap-org#1487)
Browse files Browse the repository at this point in the history
* paste Mayo's surrogate_crystallizer unit in raw form

* start code cleanup in surrogate crystallizer, add vapor prop pack config option, rename class

* git rid of self.vap_prop instances

* paste in seawater relationship for enth_mass_phase

* pass tests for enth_mass_phase implementation in MCAS

* check value of enthalpy in test

* bring in test code

* get test crystallizer surrogate file to pass

* only construct enthalpy params once

* blk

* add surrogate files

* clean up imports and other code

* track costing file

* bk

* raise config error instead of string

* remove unused imports

* changed filename and started cleanup

* make folder for crystallizer surrogate defaults to test; finish cleanup

* blk

* add folder

* cleanup imports and copyright headers

* fix file path

* updated some var names to make a bit more user friendly and closer to idaes naming; also described a config option that was unclear

* blk

* delete erronous code in mcas during initialization of enthalpy

* fix issue with file path?

* blk

* black

* move surrogate defaults

* refine test file

* add pressure_sat and add some scaling

* black

* add to MCAS tech ref; revise subtlety in sw prop model tech ref

* revise docs

* fix rst

* adjust th equation in rst

* remove units inclusion for D

* fix fraction

* fix path issue

* blk

* fix relatiosnhip

* try to fix path issues now having created init files

* revise doc

* tinker with docs

* cleanup

* add init files!

* regression tests for enth flow and pressure sat
  • Loading branch information
adam-a-a authored Oct 31, 2024
1 parent 53bd64e commit cbc8882
Show file tree
Hide file tree
Showing 14 changed files with 1,366 additions and 29 deletions.
18 changes: 14 additions & 4 deletions docs/technical_reference/property_models/mc_aq_sol.rst
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ Parameters
.. csv-table::
:header: "Description", "Symbol", "Parameter", "Index", "Units"

"Component molecular weight :sup:`1`", ":math:`m_N`", "mw_comp", "[j]", ":math:`\text{kg mol}^{-1}`"
"Component molecular weight :sup:`1`", ":math:`MW`", "mw_comp", "[j]", ":math:`\text{kg mol}^{-1}`"
"Stokes radius of solute", ":math:`r_h`", "radius_stokes_comp", "[j]", ":math:`\text{m}`"
"Molar volume of solute", ":math:`V`", "molar_volume_phase_comp", "[p, j]", ":math:`\text{m}^3 \text{ mol}^{-1}`"
"Dynamic viscosity", ":math:`\mu`", "visc_d_phase", "[p]", ":math:`\text{Pa s}`"
Expand Down Expand Up @@ -112,7 +112,11 @@ Properties
"Debye-Huckel constant A", ":math:`A`", "deby_huckel_constant", "none", ":math:`\text{dimensionless}`"
"Ionic Strength", ":math:`I`", "ionic_strength_molal", "none", ":math:`\text{mol kg}^{-1}`"
"Mass diffusivity of solute", ":math:`D`", "diffus_phase_comp", "[p, j]", ":math:`\text{m}^2 \text{ s}^{-1}`"

"Specific enthalpy", ":math:`\widehat{H}`", "enth_mass_phase", "[p]", ":math:`\text{J/kg}`"
"Enthalpy flow", ":math:`H`", "enth_flow", "None", ":math:`\text{J/s}`"
"Saturation pressure", ":math:`P_v`", "pressure_sat", "None", ":math:`\text{Pa}`"
"Total hardness as CaCO3",":math:`TH`","total_hardness", "None", ":math:`\text{mg/L}`"
"Total dissolved solids",":math:`TDS`","total_dissolved_solids", "None", ":math:`\text{mg/L}`"


Relationships
Expand All @@ -137,7 +141,13 @@ Relationships
"Phase electrical conductivity", ":math:`\lambda=\Lambda\sum_{j\in cation}{\left|z_j\right|n_j}`"
"Debye-Huckel constant", ":math:`A=\frac{\left(2 \pi N_A\right)^{0.5}}{log(10)} \left(\frac{\textbf{e}^2}{4 \pi \epsilon \epsilon_0 kT}\right)^{\frac{3}{2}}`"
"Ionic strength", ":math:`I=0.5\sum_{j\in ion}{z_j^2b_j}`"
"Component mass diffusivity :sup:`5`", ":math:`D\text{ specified in data argument}` or :math:`D \text{ }[\text{m}^2 \text{ s}^{-1}]=\frac{\chi_{1}}{(\mu \text{ }[\text{cP}])^{\chi_{2}}(V \text{ }[\text{cm}^3 \text{ mol}^{-1}])^{\chi_{3}}}`"
"Component mass diffusivity :sup:`5`", ":math:`D\text{ specified in data argument}` or :math:`D \text{ }=\frac{\chi_{1}}{(\mu \text{ }[\text{cP}])^{\chi_{2}}(V \text{ }[\text{cm}^3 \text{ mol}^{-1}])^{\chi_{3}}}`"
"Specific enthalpy", "Equations 25-27 in Nayar et al. (2016)"
"Enthalpy flow", ":math:`H = \sum_{j} M_j \cdotp \widehat{H}`"
"Saturation pressure", "Equations 5 and 6 in Nayar et al. (2016)"
"Total hardness as CaCO3",":math:`TH = \sum_{j \in \text{polyvalent cation set}} {n_j z_j} \frac{MW_{CaCO3}}{z_{CaCO3}}`"
"Total dissolved solids",":math:`TDS = \sum_{j \in \text{ion set}} m_j`"


.. note::

Expand Down Expand Up @@ -166,7 +176,7 @@ Physical/chemical constants

Scaling
-------
A comprehensive scaling factor calculation method is coded in this property package. Among the state variables (:math:`N, T, \text{and } p`), default scaling factors for :math:`T` and :math:`p` were set and do not need users' input, while, for :math:`N`, usually require a user input via an interface. The coding interface to set defalut scaling factor for :math:`N` and call the scaling calculation for other variables is the following.
A comprehensive scaling factor calculation method is coded in this property package. Among the state variables (:math:`N, T, \text{and } p`), default scaling factors for :math:`T` and :math:`p` were set and do not need users' input, while, for :math:`N`, usually require a user input via an interface. The coding interface to set default scaling factor for :math:`N` and call the scaling calculation for other variables is the following.

.. code-block::
Expand Down
12 changes: 3 additions & 9 deletions docs/technical_reference/property_models/seawater.rst
Original file line number Diff line number Diff line change
Expand Up @@ -51,18 +51,14 @@ Properties
"Latent heat of vaporization", ":math:`h_{vap}`", "dh_vap_mass", "None", ":math:`\text{J/kg}`"
"Diffusivity", ":math:`D`", "diffus_phase_comp", "[p]", ":math:`\text{m}^2\text{/s}`"
"Boiling point elevation", ":math:`BPE`", "boiling_point_elevation_phase", "[p]", ":math:`\text{K}`"


**The properties make use of the average molecular weight of sea salt, ≈ 31.40 g/mol, reported in the Reference-Composition Salinity Scale (Millero et al., 2008) to convert to moles.**

.. csv-table::
:header: "Description", "Symbol", "Variable", "Index", "Units"

"Component mole flowrate", ":math:`N_j`", "flow_mol_phase_comp", "[p, j]", ":math:`\text{mole/s}`"
"Component mole fraction", ":math:`y_j`", "mole_frac_phase_comp", "[p, j]", ":math:`\text{dimensionless}`"
"Molality", ":math:`Cm`", "molality_phase_comp", "['TDS']", ":math:`\text{mole/kg}`"
"Osmotic pressure", ":math:`\pi`", "pressure_osm_phase", "None", ":math:`\text{Pa}`"

**The properties make use of the average molecular weight of sea salt, ≈ 31.40 g/mol, reported in the Reference-Composition Salinity Scale (Millero et al., 2008) to convert to moles.**


Relationships
-------------
.. csv-table::
Expand All @@ -87,8 +83,6 @@ Relationships
"Diffusivity", "Equation 6 in Bartholomew et al. (2019)"
"Boiling point elevation", "Equation 36 in Sharqawy et al. (2010)"



Note: Osmotic pressure calculation (based on equation 48 in Nayar et al. (2016)) uses the density of water as a function of temperature (:math:`\rho_w`) and the ideal gas constant (:math:`R\text{, 8.314 J/mol}\cdotp\text{K}`), in addition to previously defined variables.

Scaling
Expand Down
189 changes: 189 additions & 0 deletions watertap/costing/unit_models/surrogate_crystallizer.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,189 @@
#################################################################################
# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California,
# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory,
# National Renewable Energy Laboratory, and National Energy Technology
# Laboratory (subject to receipt of any required approvals from the U.S. Dept.
# of Energy). All rights reserved.
#
# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license
# information, respectively. These files are also available online at the URL
# "https://github.com/watertap-org/watertap/"
#################################################################################

import pyomo.environ as pyo
from watertap.costing.util import (
register_costing_parameter_block,
make_capital_cost_var,
)


def build_surrogate_crystallizer_cost_param_block(blk):

blk.steam_pressure = pyo.Var(
initialize=3,
units=pyo.units.bar,
doc="Steam pressure (gauge) for crystallizer heating: 3 bar default based on Dutta example",
)

blk.efficiency_pump = pyo.Var(
initialize=0.7,
units=pyo.units.dimensionless,
doc="Crystallizer pump efficiency - assumed",
)

blk.pump_head_height = pyo.Var(
initialize=1,
units=pyo.units.m,
doc="Crystallizer pump head height - assumed, unvalidated",
)

# Crystallizer operating cost information from literature
blk.fob_unit_cost = pyo.Var(
initialize=675000,
doc="Forced circulation crystallizer reference free-on-board cost (Woods, 2007)",
units=pyo.units.USD_2007,
)

blk.ref_capacity = pyo.Var(
initialize=1,
doc="Forced circulation crystallizer reference crystal capacity (Woods, 2007)",
units=pyo.units.kg / pyo.units.s,
)

blk.ref_exponent = pyo.Var(
initialize=0.53,
doc="Forced circulation crystallizer cost exponent factor (Woods, 2007)",
units=pyo.units.dimensionless,
)

blk.iec_percent = pyo.Var(
initialize=1.43,
doc="Forced circulation crystallizer installed equipment cost (Diab and Gerogiorgis, 2017)",
units=pyo.units.dimensionless,
)

blk.steam_cost = pyo.Var(
initialize=0.004,
units=pyo.units.USD_2018 / (pyo.units.meter**3),
doc="Steam cost, Panagopoulos (2019)",
)

costing = blk.parent_block()
costing.register_flow_type("steam", blk.steam_cost)


def cost_surrogate_crystallizer(blk):
"""
Function for costing the surrogate crystallizer by the mass flow of produced crystals.
The operating cost model assumes that heat is supplied via condensation of saturated steam (see Dutta et al.)
"""
cost_crystallizer_by_crystal_mass(blk)


def _cost_crystallizer_flows(blk):
blk.costing_package.cost_flow(
pyo.units.convert(
(blk.unit_model.heat_required / _compute_steam_properties(blk)),
to_units=pyo.units.m**3 / pyo.units.s,
),
"steam",
)


@register_costing_parameter_block(
build_rule=build_surrogate_crystallizer_cost_param_block,
parameter_block_name="surrogate_crystallizer",
)
def cost_crystallizer_by_crystal_mass(blk):
"""
Mass-based capital cost for FC crystallizer
"""
make_capital_cost_var(blk)
blk.costing_package.add_cost_factor(blk, "TIC")
blk.cost_factor = 1 # blk.costing_package.add_cost_factor(blk, "TIC")
blk.capital_cost_constraint = pyo.Constraint(
expr=blk.capital_cost
== blk.cost_factor
* pyo.units.convert(
(
blk.costing_package.surrogate_crystallizer.iec_percent
* blk.costing_package.surrogate_crystallizer.fob_unit_cost
* (
blk.unit_model.flow_mass_sol_total
/ blk.costing_package.surrogate_crystallizer.ref_capacity
)
** blk.costing_package.surrogate_crystallizer.ref_exponent
),
to_units=blk.costing_package.base_currency,
)
)
_cost_crystallizer_flows(blk)


def _compute_steam_properties(blk):
"""
Function for computing saturated steam properties for thermal heating estimation.
Args:
pressure_sat: Steam gauge pressure in bar
Out:
Steam thermal capacity (latent heat of condensation * density) in kJ/m3
"""
pressure_sat = blk.costing_package.surrogate_crystallizer.steam_pressure
# 1. Compute saturation temperature of steam: computed from El-Dessouky expression
tsat_constants = [
42.6776 * pyo.units.K,
-3892.7 * pyo.units.K,
1000 * pyo.units.kPa,
-9.48654 * pyo.units.dimensionless,
]
psat = (
pyo.units.convert(pressure_sat, to_units=pyo.units.kPa)
+ 101.325 * pyo.units.kPa
)
temperature_sat = tsat_constants[0] + tsat_constants[1] / (
pyo.log(psat / tsat_constants[2]) + tsat_constants[3]
)

# 2. Compute latent heat of condensation/vaporization: computed from Sharqawy expression
t = temperature_sat - 273.15 * pyo.units.K
enth_mass_units = pyo.units.J / pyo.units.kg
t_inv_units = pyo.units.K**-1
dh_constants = [
2.501e6 * enth_mass_units,
-2.369e3 * enth_mass_units * t_inv_units**1,
2.678e-1 * enth_mass_units * t_inv_units**2,
-8.103e-3 * enth_mass_units * t_inv_units**3,
-2.079e-5 * enth_mass_units * t_inv_units**4,
]
dh_vap = (
dh_constants[0]
+ dh_constants[1] * t
+ dh_constants[2] * t**2
+ dh_constants[3] * t**3
+ dh_constants[4] * t**4
)
dh_vap = pyo.units.convert(dh_vap, to_units=pyo.units.kJ / pyo.units.kg)

# 3. Compute specific volume: computed from Affandi expression (Eq 5)
t_critical = 647.096 * pyo.units.K
t_red = temperature_sat / t_critical # Reduced temperature
sp_vol_constants = [
-7.75883 * pyo.units.dimensionless,
3.23753 * pyo.units.dimensionless,
2.05755 * pyo.units.dimensionless,
-0.06052 * pyo.units.dimensionless,
0.00529 * pyo.units.dimensionless,
]
log_sp_vol = (
sp_vol_constants[0]
+ sp_vol_constants[1] * (pyo.log(1 / t_red)) ** 0.4
+ sp_vol_constants[2] / (t_red**2)
+ sp_vol_constants[3] / (t_red**4)
+ sp_vol_constants[4] / (t_red**5)
)
sp_vol = pyo.exp(log_sp_vol) * pyo.units.m**3 / pyo.units.kg

# 4. Return specific energy: density * latent heat
return dh_vap / sp_vol
Empty file.
Loading

0 comments on commit cbc8882

Please sign in to comment.